Re: [Gretl-users] Kalman Filter
On Thu, March 16, 2006 13:29, Ignacio Díaz-Emparanza wrote:
Attached you will find a script that applies a Hodrick-Prescott
filter by
means of a Kalman filter. I hope in future versions of gretl the "function"
command can accept matrices as parameters, so this will do easier to apply
the KF. We could advance then in programming, for example, Maximum Likelihood
trough "prediction error decomposition" and then Structural Time Series
Models.
With current CVS you already can have matrices as arguments. What you can't do
(and I admit that needs to be addressed) is *return* matrices from functions.
The following is an example function, basically ripped from your script with
some optimisations, which does the forward pass and returns the state as a
series.
function kfilt(series y, matrix a0, matrix p0, matrix bT, matrix Z, \
matrix bT, matrix Sigma_w, scalar sigma_e)
# Forward solution
#
scalar T = rows(Z)
scalar m = rows(a0)
matrix at_t = a0
matrix pt_t = p0
matrix at = zeros(m,T)
printf "\nStarting the filter process\n"
printf "\n...Filtering...\n"
loop for i=1..T
# Prediction equations
# eq. (2.2a)
matrix at_t = bT*at_t
# eq. (2.2b)
matrix pt_t1 = bT*pt_t*transp(bT)+Sigma_w
# eq. (2.3c)
matrix zt = Z[$i,]
matrix H = pt_t1*zt'
matrix f = zt*H + sigma_e
if i>1
matrix pstar_t=pt_t*transp(bT)*inv(pt_t1)
endif
#Updating equations
# eq. (2.4a)
scalar e = y[$i] - zt*at_t
matrix at_t = at_t + H*(e/f)
# eq (2.4b)
matrix pt_t = pt_t1 - H*H' * (1/f)
matrix at[,$i]=at_t
if i==1
matrix pt=pt_t
else
if i==2
matrix pt=pt~pt_t
matrix pstar=pstar_t
else
matrix pt=pt~pt_t
matrix pstar=pstar~pstar_t
endif
endif
end loop
series filtered = at[1,]
printf "\nFilter done\n"
return series filtered
end function
Riccardo "Jack" Lucchetti
Dipartimento di Economia
Facoltà di Economia "G. Fuà"
Ancona