You have done a great job sir. How I wish I could contribute to this project but my knowledge in hansl is very limited. I have a dream of developing a GUI package to estimate DSGE models in gretl due to it user friendly nature. Please, if possible that you later improve this code, do not hesitate to post it here. 

Good luck.

On Sun, Feb 23, 2020, 3:30 PM Mario Marchetti <marchetti224@gmail.com> wrote:
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@gretlml.univpm.it
To unsubscribe send an email to gretl-users-leave@gretlml.univpm.it
Website: https://gretlml.univpm.it/postorius/lists/gretl-users.gretlml.univpm.it/