A simple Real Business Cycle model with gretl
by Mario Marchetti
Goodmorning everyone,
I'm Mario and for fun and study, especially to practice in hansl, I am trying to adapt a Maltlab script to the hansl language.
This script was written by Ryo Kato and consists of solving a simple RBC model.
The code is available here: http://www.ryokato.org/genmac/RBC1.m
A first (spartan) draft of the code that I wrote in hansl is the following (also available on github: https://github.com/mariometrics/RBCgretl):
<hansl>
####------------------------------------------------------------------------#####
set echo off
set messages off
## Mario Marchetti 23-02-2020
## Basic RBC model ##
## Adapted in hansl language from the code written in Matlab by Ryo Kato in 2004
## ------------------- [1] Parameter proc ------------------------
sigma = 1.5 # CRRA
alpha = 0.3 # Cobb-Dag
myu = 1 # labor-consumption supply
beta = 0.99 # discount factor
delta = 0.025 #depreciation
lamda = 2 # labor supply elasticity >1
phi = 0.8 # AR(1) in tech
param = {sigma,alpha,myu,beta,delta,lamda,phi}
## --------------------- [2] Steady State proc >> -----------------------
# SS capital & ss labor
# (1) real rate (By SS euler)
kls = (((1/beta)+delta-1)/alpha)^(1/(alpha-1))
# (2) wage
wstar = (1-alpha)*(kls)^alpha
# (3) Labor and goods market clear
clstar = kls^alpha - delta*kls
lstar = ((wstar/myu)*(clstar^(-sigma)))^(1/(lamda+sigma))
kstar = kls*lstar
cstar = clstar*lstar
vstar = 1
Ystar = (kstar^alpha)*(lstar^(1-alpha))
ssCKoLY = {cstar,kstar;lstar,Ystar} # show SS values
## --------------------------[2] MODEL proc-----------------------------##
function matrix RBC(matrix *param,matrix *x)
sigma = param[1]
alpha = param[2]
myu = param[3]
beta = param[4]
delta = param[5]
lamda = param[6]
phi = param[7]
# Define endogenous vars ('a' denotes t+1 values)
la = x[1]
ca = x[2]
ka = x[3]
va = x[4]
lt = x[5]
ct = x[6]
kt = x[7]
vt = x[8]
ra = 0
rt = 0
# Eliminate Price
ra = (va*alpha*(ka/la)^(alpha-1))
wt = (1-alpha)*vt*(kt/lt)^alpha
# Optimal Conditions & state transition
labor = lt^lamda-wt/(myu*ct^sigma) # LS = LD
euler = ct^(-sigma) -(ca^(-sigma))*beta*(1+ra-delta) # C-Euler
capital = ka - (1-delta)*kt-vt*(kt^alpha)*(lt^(1-alpha))+ct # K-trans
tech = va - phi*vt
matrix optcon = {labor;euler;capital;tech}
return optcon
end function
function scalar RBCY(matrix *param,matrix *xr)
# GDP (Optional)
alpha = param[2]
vt = xr[3]
kt = xr[2]
lt = xr[1]
Yt = vt*(kt^alpha)*(lt^(1-alpha))
return Yt
end function
# Evaluate each derivate
matrix x = {lstar,cstar,kstar,vstar,lstar,cstar,kstar,vstar}
matrix xr = {lstar,kstar,vstar}
# Numerical jacobian
matrix coeff = fdjac(x,RBC(¶m,&x))
matrix coeffy = fdjac(xr,RBCY(¶m,&xr))
# In terms of # deviations from ss
matrix vo = {lstar,cstar,kstar,vstar}
matrix TW = vo | vo | vo | vo
matrix B = -coeff[,1:4].*TW
matrix C = coeff[,5:8].*TW
# B[c(t+1) l(t+1) k(t+1) z(t+1)] = C[c(t) l(t) k(t) z(t)]
matrix A = inv(C)*B #(Linearized reduced form )
# For GDP( optional)
matrix ve = {lstar,kstar,vstar}
matrix NOM = {Ystar,Ystar,Ystar}
matrix PPX = coeffy.*ve./NOM
## =========== [4] Solution proc ============== ##
# EIGEN DECOMPOSITION
matrix W = {}
matrix theta = eigengen(A, &W)
Q = inv(W)
V = zeros(4,4)
V[diag] = theta
LL = W*V*Q # not find a role yet...
# Extract stable vectors
matrix SQ = {}
loop j = 1..rows(theta) --quiet
if abs(theta[j]) > 1.000000001
SQ |= Q[j,]
endif
endloop
# Extract unstable vectors
matrix UQ = {}
loop jj = 1..rows(theta) --quiet
if abs(theta[jj])<0.9999999999
UQ |= Q[jj,]
endif
endloop
# Extract stable roots
matrix VLL = {}
loop jjj = 1..rows(theta) --quiet
if abs(theta[jjj]) >1.0000000001
VLL |= theta[jjj,]
endif
endloop
# [3] ELIMINATING UNSTABLE VECTORS
k = min({rows(SQ),cols(SQ)}) # # of predetermined vars
n = min({rows(UQ),cols(UQ)}) # # of jump vars
nk = {n,k}
# Stable V (eig mat)
diago = zeros(rows(VLL),rows(VLL))
diago[diag] = VLL
VL = inv(diago)
# Elements in Q
PA = UQ[1:n,1:n]
PB = UQ[1:n,n+1:n+k]
PC = SQ[1:k,1:n]
PD = SQ[1:k,n+1:n+k]
P = -inv(PA)*PB # X(t) = P*S(t)
PE = PC*P+PD
# SOLUTION
PX = inv(PE)*VL*PE
AA = Re(PX)
## ------------------ [5] SIMULATION proc ----------------- ##
# [4] TIME&INITIAL VALUES
t = 48 # Time span
# Initial Values
# state var + e
S1 = {0;0.06}
# [5] SIMULATION
Ss = S1
S = zeros(t,k)
loop i = 1..t --quiet
q = AA*Ss
S[i,] = q'
Ss = S[i,]'
endloop
SY = S1' | S
X = (Re(P)*SY')'
# Re-definition
ci = X[,1]
li = X[,2]
ki = SY[,1]
vi = SY[,2]
Yi = (PPX*XI')'
# [6] DRAWING FIGURES
gnuplot --matrix=Yi --time-series --with-lines --output=display { set linetype 3 lc rgb "#0000ff"; set title "Y"; set key rmargin; set xlabel "time"; set ylabel "IRF Y_t"; }
# put columns together and add labels
plotmat = X ~ SY
strings cnames = defarray("C", "L","K","V")
cnameset(plotmat, cnames)
scatters 1 2 3 4 --matrix=plotmat --with-lines --output=display
####-------------------------------------------------------------------#####
</hansl>
So I write to ask for suggestions to improve or "streamline" the code, or to help find errors that have escaped me.
All to try to improve my knowledge of gretl software and its scripting language: "hansl".
For example: how can I improve the jacobian calculation?
Thanks to everyone and have a good day.
6 years, 3 months
Reading error from excel files created by Ox
by Fred Engst
Hi Allin and all other hard working member of the gretl team.
I’m having a hard time reading excel files created by "Ox Console version 8.02 (OS_X_64/U) (C) J.A. Doornik, 1994-2018”.
If I save a matrix in xls format, gretl gives me a message of “Failed to get workbook info”.
If I save the matrix in xlsx format, gretl skips the header and gives me generic variable names as in: v1, v2, …
Any suggestion for what I should do?
Fred
6 years, 3 months
MLE and binary (scalar) min and max operators
by Alecos Papadopoulos
Good evening. Will the mle command in gretl have any compatibility
problem if in the likelihood some of the parameters under estimation
appear also inside binary min and max operators
Namely, a thing like (bogus likelihood),
<hansl>
scalar a = starting value
scalar b = starting value
scalar s = starting value
mle logl = log(a) + log(Φ(ε / s + min(a-b,0))) + exp(-max(a,b))
...
</hansl>
I am not asking about convergence issues or negative logarithms, these
are model/data issues. Only if the mle command accepts in principle to
work with the binary min/max operators.
--
Alecos Papadopoulos PhD
Athens University of Economics and Business
web: alecospapadopoulos.wordpress.com/
skype:alecos.papadopoulos
6 years, 3 months
Identification of SVARs using heteroskedasticity
by anzervas＠yahoo.com
Dear all (especially Sven and Riccardo),
A recent strand in SVAR literature uses heteroskedasticity to identify the structural shocks - those interested in the topic may read mainly Rigobon (Review of Economics and Statistics 2003), Lanne and Lutkepohl (Journal of Money Credit and Banking 2008) and / or Bacchiocchi and Fanelli (Oxford Bulletin of Economics and Statistics 2015). If one has 2 variance - covariance matrices, there is no need to add zero or other restrictions to identify the structural matrices A or B.
It is not clear to me how one can implement this on Gretl. Especially I did not find any way to implement this method using matrix operations, and the only possible way seems to employ an optimization algorithm. I have tried it, but in vain. Below I attach a script, formulated like in Lanne - Lutkepohl (2008) but it is not working properly. Especially, I cannot understand how to make the procedure go to a solution that respects the equality S1 = B*B'.
In addition, Rigobon (op. cit.) mentions that he uses GMM to solve the problem. I do not know if this is feasible in Gretl, but even if it is, it is not clear to me how one could write the GMM block to do it (though the manual mentions that one may use only matrices in GMM block equations).
Any suggestions / corrections are welcome. This also seems be a good functionality to add to the SVAR package (in addition to add AB model functionality to VECMs, for completeness purposes).
Kind regards,
Andreas
<begin hansl scipt>
set verbose off
# set lbfgs on
set bfgs_maxiter 50000
set bfgs_toler 0.0000000000000001
function scalar LnL(const matrix param, matrix S1, matrix S2)
A = {param[1],param[2],param[3],param[4]; \
param[5],param[6],param[6],param[8]; \
param[9],param[10],param[11],param[12]; \
param[13],param[14],param[15],param[16]}
L = zeros(rows(S1),cols(S1))
L[1,1] = param[17]
L[2,2] = param[18]
L[3,3] = param[19]
L[4,4] = param[20]
LL0 = -100*0.5*(ln(det(S1)) + tr(S1*inv(S1))) \
-100*0.5*(ln(det(S2)) + tr(S2*inv(S2)))
LLi = -100*0.5*(ln(det(A*A')) + tr(S1*inv(A*A'))) \
-100*0.5*(ln(det(A*L*A')) + tr(S2*inv(A*L*A')))
# dist = abs(LLi - LL0)
dist = maxc(abs((vech(S1)|vech(S2)) - (vech(A*A')|vech(A*L*A'))))
dist
return dist
end function
# structural shocks and matrix
E1 = mnormal(100,4)
E2 = mnormal(100,cols(E1)).*{0.1, 2.5, 0.4, 1.7}
W = 0.1*mrandgen(i, 1, 4, 4, 4)
# reduced form residuals
U1 = E1*W
U2 = E2*W
S1 = U1'U1/rows(U1)
S2 = U2'U2/rows(U2)
# initial parameters
param = 0.1*abs(mnormal(cols(U1)^2+cols(U1),1))
# minimization
ff=BFGSmin(¶m, LnL(param, S1, S2))
# bounds = seq(1,rows(param))'~ones(rows(param),2).*{-20, 20}
# bounds[17,] = {17, 0, 10}
# bounds[18,] = {18, 0, 10}
# bounds[19,] = {19, 0, 10}
# bounds[20,] = {20, 0, 10}
# ffc=BFGScmin(¶m, bounds, LnL(param, S1, S2))
param
B = {param[1],param[2],param[3],param[4]; \
param[5],param[6],param[6],param[8]; \
param[9],param[10],param[11],param[12]; \
param[13],param[14],param[15],param[16]}
L = zeros(rows(S2),cols(S2))
L[1,1] = param[17]
L[2,2] = param[18]
L[3,3] = param[19]
L[4,4] = param[20]
S1b = B*B'
S2b = B*L*B'
# Check that estimations reproduce reduced form covariance matrices and structural matrix
S1
S1b
L
S2
S2b
W
B
<end hansl scipt>
6 years, 3 months
random-effects Tobit model
by Sven Schreiber
Hi,
almost two years ago it was asked on this mailing list how to estimate a
Tobit model on panel data with gretl.
It has taken a while, but now there's a function package "retobit" to do
just that. As with all (moderated) function packages this is installable
directly from within gretl from the package server.
You can see the help text before installing by going to
http://ricardo.ecn.wfu.edu/gretl/cgi-bin/gretldata.cgi?opt=SHOW_FUNCS
and then clicking on "retobit.gfn" (currently number 98 in the
alphabetical list).
Feedback (including problem reports) is very welcome.
cheers
sven
6 years, 3 months
Gretl from USB stick under windows
by Stefano Fachin
I can confirm that (as expected) Gretl works smoothly from a USB stick.
I had a class of 50 students using it on Linux PCs running a Windows
virtual machine, no problems
bye
Stefano
--
_________________________________________________________________________
Stefano Fachin
Professore Ordinario di Statistica Economica
Dip. di Scienze Statistiche
Università di Roma "La Sapienza"
P.le A. Moro 5 - 00185 Roma - Italia
Tel. +39-06-49910834
fax +39-06-49910072
web http://stefanofachin.site.uniroma1.it/
--
________________________________________________________
Le informazioni
contenute in questo messaggio di posta elettronica sono strettamente
riservate e indirizzate esclusivamente al destinatario. Si prega di non
leggere, fare copia, inoltrare a terzi o conservare tale messaggio se non
si è il legittimo destinatario dello stesso. Qualora tale messaggio sia
stato ricevuto per errore, si prega di restituirlo al mittente e di
cancellarlo permanentemente dal proprio computer.
The information
contained in this e mail message is strictly confidential and intended for
the use of the addressee only. If you are not the intended recipient,
please do not read, copy, forward or store it on your computer. If you have
received the message in error, please forward it back to the sender and
delete it permanently from your computer system.
6 years, 3 months
Mac woes
by Riccardo (Jack) Lucchetti
Sorry folks,
back in December, a few users wrote me about problems that they were
having when opening files on their Macs; one of them sent me the attached
screenshot. Apparently, the gretl app has no access to the Desktop.
At the time, I was busy with other stuff but then it completely slipped my
mind. I'm no mac user, so I really can't say what's going on and if the
problem is still there. I tried to search the list for similar problems
but found nothing. Any hints?
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
6 years, 3 months
heckprobit package updated
by Sven Schreiber
Hello microeconometricians out there,
the contributed function package heckprobit.gfn by Claudia Pigini has
been updated to version 0.2, which proves that even packages that
haven't been updated for years don't have to be dead necessarily.
(As the name suggests, the package offers estimation of Probit models
with Heckman-style sample selection.)
The changes are mostly internal, but they now require a gretl version of
2017b or later. So if you're a user of that package and haven't upgraded
your gretl installation in years, now's the time.
cheers
sven
6 years, 3 months
dbnomics bundle returns empty list
by klaus.hasenbach＠web.de
Hallo,
not with all but in some databases my old bundle requests does not work
anymore.
For instance:
open dbnomics
set verbose on
include dbnomics.gfn
nulldata 900
setobs 12 1948.1
# OECD Deflators: "countries.typeGDP.DOBSA.(Q)uarterly"
provider = "OECD"
database = "QNA"
bundle spec = defbundle("mask","G-20+DEU+USA+FRA+JPN+CHN+CHE+OECD..DOBSA.Q")
bs = dbnomics_get_multiple(provider, database, 1000, 0, spec)
dbnomics_bundles_print(bs)
list X = dbnomics_bundles_to_list( bs, "series_code" )
printf "\nHere are the series in list X:\n"
list X print
or
# Bis Dept:
"(Q)uartely.countries.debtors.creditors.(M)arketvalue.(770)%GDP.A"
provider = "BIS"
database = "total_credit"
bundle spec = defbundle("mask","Q....M.770.A")
bs = dbnomics_get_multiple(provider, database, 1000, 0, spec)
dbnomics_bundles_print(bs)
list X = dbnomics_bundles_to_list( bs, "series_code" )
printf "\nHere are the series in list X:\n"
list X print
Gretl finds the series but it returns "empty list" for the observations
, using either dbnomics v0.2 or v0.35 and running gretl with MS Windows.
(Fetching those series with menue option in gretl works fine.)
Is there a way to solve this? Thank you.
Klaus Hasenbach
6 years, 3 months
