[Gretl-users] Re: A simple Real Business Cycle model with gretl

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(&param,&x)) matrix coeffy = fdjac(xr,RBCY(&param,&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. _______________________________________________ Gretl-users mailing list -- gretl-users(a)gretlml.univpm.it To unsubscribe send an email to gretl-users-leave(a)gretlml.univpm.it Website: https://gretlml.univpm.it/postorius/lists/gretl-users.gretlml.univpm.it/