[Gretl-devel] Problem with a package
Ignacio Diaz-Emparanza
ignacio.diaz-emparanza at ehu.es
Tue Jun 5 06:20:22 EDT 2007
I am building some packages to estimate Structural Time Series Models. Now I
have a collection of functions that work well in the console but I am having
problems in packaging them. An example is the package enclosed.
I can run it in the console without any problem, for example:
(I know it is incorrect to use a series like this because of the seasonality,
but it is very illustrative)
-------
data9-3.gdt
include LLTestim.gfn
list comp = LLTestim(reskwh)
-------
(After some time...) This gives a table output and a "comp" list with two
series LLTtrend, and LLTslope. But I cannot do it to work on the GUI,
I always get an "Unespecified error". Any idea on how to solve this?
--
Ignacio Diaz-Emparanza
DEPARTAMENTO DE ECONOMÍA APLICADA III (ECONOMETRÍA Y ESTADÍSTICA)
UPV/EHU
Avda. Lehendakari Aguirre, 83 | 48015 BILBAO
T.: +34 946013732 | F.: +34 946013754
www.et.bs.ehu.es
-------------- next part --------------
<?xml version="1.0" encoding="UTF-8"?>
<gretl-functions>
<gretl-function-package name="LLTestim" ID="1181037704" needs-time-series-data="true" minver="1.6.5">
<author>Ignacio Diaz-Emparanza</author>
<version>1.0</version>
<date>2007-06-05</date>
<description>Local Linear Trend Model</description>
<gretl-function name="LLTestim" private="0">
<params count="6">
<param name="y" type="series"/>
<param name="q" type="matrixref" optional="true"/>
<param name="irreg" type="bool" default="1"/>
<param name="level" type="bool" default="1"/>
<param name="slope" type="bool" default="1"/>
<param name="sigmatol" type="scalar" default="0.0001"/>
</params>
<return name="compo" type="list"/>
<code>if (irreg!=1&&irreg!=0) || (level!=1&&level!=0) || (slope!=1&&slope!=0)
funcerr "irreg, level and slope must be 0 or 1"
end if
if $pd>1
print "LLTestim warning: your data are seasonal, you should use the BSMestim function"
end if
matrix fijos = { irreg, level, slope }
matrix theta = { -0.5, -1.5}
matrix at
series V = 0
series F = 0
scalar logLc=0
matrix pstar
scalar sigma_e=1
matrix bT
set bfgs_toler 1.E-2
M = BFGSmax(theta, "LLT(&theta, y, sigma_e, &bT, &at, &pstar, &V, &F, &logLc, sigmatol, 1, fijos)")
scalar sstart = int(min(t))
scalar send = int(max(t))
scalar T=send-sstart+1
scalar sstar = (1/T)*(sum((V^2)/F)
scalar q1=exp(2*theta[1])*sstar
scalar q2=exp(2*theta[2])*sstar
matrix qp = { sstar, q1, q2 }
matrix conc=imaxr(qp)
scalar concent = conc[1]
set bfgs_toler default
matrix theta = { -0.5, -1.5 }
M = BFGSmax(theta, "LLT(&theta, y, sigma_e, &bT, &at, &pstar, &V, &F, &logLc, sigmatol, concent, fijos)")
scalar sstar = (1/T)*(sum((V^2)/F)
scalar q1=exp(2*theta[1])
scalar q2=exp(2*theta[2])
series LLTtrend = kf_smooth(pstar, &at, bT)
series LLTslope = at[2,]
list compo = LLTtrend LLTslope
printf "\nLocal Linear Trend Model estimation:\n"
printf "-----------------------------------------\n"
if concent=1
printf " sigma*=\t\t Var(eps)=%8.6E\n", sstar
printf " q1=%8.5f,\t Var(eta)=%8.6E\n", q1, q1*sstar
printf " q2=%8.5f,\t Var(xi)=%8.6E\n", q2, q2*sstar
matrix q = {1, q1, q2}'
else
if concent=2
printf " q1=%8.5f,\t Var(eps)=%8.6E\n", q1, q1*sstar
printf " sigma*=\t\t Var(eta)=%8.6E\n", sstar
printf " q2=%8.5f,\t Var(xi)=%8.6E\n", q2, q2*sstar
matrix q = {q1, 1, q2}'
else
if concent=3
printf "q1=%8.5f,\t Var(eps)=%8.6E\n", q1, q1*sstar
printf "q2=%8.5f,\t Var(eta)=%8.6E\n", q2, q2*sstar
printf "sigma*=\t\t Var(xi)=%8.6E\n", sstar
matrix q = {q1, q2, 1}'
end if
end if
end if
printf "------------------------------------------\n \n"
</code>
</gretl-function>
<gretl-function name="kf_filt" private="1">
<params count="12">
<param name="y" type="series"/>
<param name="a0" type="matrix"/>
<param name="p0" type="matrix"/>
<param name="bT" type="matrix"/>
<param name="Z" type="matrix"/>
<param name="Sigma_w" type="matrix"/>
<param name="sigma_e" type="scalar"/>
<param name="at" type="matrixref"/>
<param name="pstar" type="matrixref"/>
<param name="V" type="seriesref"/>
<param name="F" type="seriesref"/>
<param name="logLc" type="scalarref"/>
</params>
<return name="filtered" type="series"/>
<code>/*
Measurement equation:
y(t) = Z[t,]*alpha(t)+e(t) (1.1a)
State transition:
alpha(t)=bT*alpha(t-1)+w(t); (1.2a)
bT is for "bold T" and w(t)=R(t)*eta(t) in Harvey 1990
a(t) is the estimator of alpha(t)
Parameters:
y = observable series
a0 = m x 1 vector, prior a(0)
p0 = m x m matrix, prior p(0)=var(a(0))
bT = m x m matrix (transition matrix)
Z = T x m matrix
Sigma_w = m x m symmetric matrix of variance of w(t), fixed for all t
sigma_e = scalar variance of e(t), fixed for all t
at = m x T matrix (output) with the estimated states
pstar = m^2 x T matrix (output)
*/
#
# Forward solution
#
scalar T = rows(Z)
scalar m = rows(a0)
matrix at_t = a0
matrix pt_t = p0
matrix at = zeros(m,T)
# printf "\n...Filtering...\n"
loop for i=1..T --quiet
# Prediction equations
# eq. (2.2a)
matrix at_t = bT*at_t
# eq. (2.2b)
matrix pt_t1 = qform(bT,pt_t)+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 = (bT*pt_t)' inv(pt_t1)
endif
#Updating equations
# eq. (2.4a)
genr V[$i] = y[$i] - zt*at_t
matrix at_t = at_t + H*(V[$i]/f)
genr F[$i] = f
# eq (2.4b)
matrix pt_t = pt_t1 - H*H' * (1/f)
matrix at[,$i]=at_t
if i>1
if i=2
matrix pstar=vec(pstar_t)
else
matrix pstar=pstar~vec(pstar_t)
endif
endif
end loop
#Concentrated Log-likelihood fuction
scalar logLc = -(T/2)*(log(2*pi)+1)-(1/2)*sum(log(F))-(T/2)*log((1/T)*sum((V^2)/F)
series filtered = at[1,]
# printf "\nFilter done\n"
</code>
</gretl-function>
<gretl-function name="kf_smooth" private="1">
<params count="3">
<param name="pstar" type="matrix"/>
<param name="at" type="matrixref"/>
<param name="bT" type="matrix"/>
</params>
<return name="ret" type="series"/>
<code>/*
Fixed-interval smoothing
The matrix pt is not used here, but could be used if
one wants to calculate confidence intervals for the
at estimators.
*/
scalar m = rows(at)
scalar T = cols(at)
#printf "\n...Smoothing...\n"
scalar T1=T-1
loop for i=1..T1 --quiet
scalar j=T-i
matrix pstar_t = mshape(pstar[,j],m,m)
# eq. (2.9a)
matrix at[,j] += pstar_t*(at[,(j+1)]-bT*at[,j])
end loop
series ret = at[1,]
#printf "\nSmoothing done\n"
</code>
</gretl-function>
<gretl-function name="LLT" private="1">
<params count="12">
<param name="param" type="matrixref"/>
<param name="y" type="series"/>
<param name="sigma_e" type="scalar"/>
<param name="bT" type="matrixref"/>
<param name="at" type="matrixref"/>
<param name="pstar" type="matrixref"/>
<param name="V" type="seriesref"/>
<param name="F" type="seriesref"/>
<param name="logLc" type="scalarref"/>
<param name="sigmatol" type="scalar" default="0.0001"/>
<param name="concent" type="scalar" default="1"/>
<param name="fixed" type="matrix"/>
</params>
<return name="logLc" type="scalar"/>
<code>if concent>4
funcerr "concent must be 1, 2, or 3"
end if
genr time
scalar sstart = int(min(time))
scalar send = int(max(time))
scalar T=send-sstart+1
matrix bT = { 1, 1; 0, 1 }
scalar m = cols(bT)
matrix Z = ones(T,1) ~ zeros(T,1)
matrix a0 = y[sstart] * ones(2,1)
matrix p0 = 400000*I(2)
matrix Sigma_w = zeros(2,2)
if concent=1
scalar sigma_e=1
scalar tmp = exp(2*param[1])*fixed[2]
Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
matrix param[1] = (tmp>sigmatol) ? param[1] : -500
scalar tmp = exp(2*param[2])*fixed[3]
Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
matrix param[2] = (tmp>sigmatol) ? param[2] : -500
else
if concent=2
scalar tmp = exp(2*param[1])*fixed[1]
scalar sigma_e=(tmp>sigmatol) ? tmp : 0
matrix param[1] = (tmp>sigmatol) ? param[1] : -500
Sigma_w[1,1] = 1
scalar tmp = exp(2*param[2])*fixed[3]
Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
matrix param[2] = (tmp>sigmatol) ? param[2] : -500
else
if concent=3
scalar tmp = exp(2*param[1])*fixed[1]
scalar sigma_e=(tmp>sigmatol) ? tmp : 0
matrix param[1] = (tmp>sigmatol) ? param[1] : -500
Sigma_w[2,2] = 1
scalar tmp = exp(2*param[2])*fixed[2]
Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
matrix param[2] = (tmp>sigmatol) ? param[2] : -500
endif
endif
endif
series filtered = kf_filt(y, a0, p0, bT, Z, Sigma_w, sigma_e, &at, &pstar, &V, &F, &logLc)
</code>
</gretl-function>
</gretl-function-package>
</gretl-functions>
More information about the Gretl-devel
mailing list