1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
|
/*
BSMestim.inp
The code in this script comes from the retired function package "BSMestim"
by Ignacio Diaz-Emparanza, which has been superseded by the "StrucTiSM"
package.
It shows how a Kalman filter can be programmed "by hand" in gretl's
language hansl. Nowadays this is not necessary anymore, since gretl has
built-in Kalman functionality, but Ignacio's script may still be
interesting for pedagogical reasons.
Below this comment the original help text of the package is reproduced,
then up to about line 365 the necessary functions are defined, then the
code from the old package's sample script appears below.
(Some syntax has been marginally edited by Sven.)
*/
/*
Help text:
This package estimates a Basic Structural Model as defined in Harvey (1991)
"Time Series Models".
The model is writen in state-space form, and the program uses the prediction
error decomposition form of the likelihood, and evaluates the likelihood
function using the Kalman Filter. The gretl function "BFGSmax" is used to
obtain the parameters that maximize the likelihood function.
Basic Structural Model
----------------------
The model:
y(t) = mu(t) + gamma(t) + e(t)
mu(t) = mu(t-1)+beta(t-1)+eta(t)
beta(t) = beta(t-1)+xi(t)
gamma(t) = -SUM_j=1^s-1 gamma(t-j) + omega(t)
in state-space form:
y(t) = Z(t)*alpha(t) + e(t)
alpha(t) = bT*alpha(t-1) + w(t)
where
alpha(t) = [mu(t), beta(t), gamma(t), ..., gamma(t-(s-2))]'
bT = ( 1, 1, 0, ... , 0
0, 1, 0, ... , 0
0, 0, -1, -1, .., -1
0, 0, 1, 0, ..., 0
.............
0, 0, 0, 0,.. 1, 0)
Z(t) = (1, 0, 1, 0, ... ,0) for all t
w(t) = [eta(t), xi(t), omega(t), 0, ..., 0]'
Sigma_w = sigma^2_e (q1, 0, .. 0
0, q2, ..0
0, 0, q3, .. 0
0, 0, 0, ... 0
.............
0, 0, 0, ... 0 )
--------------------------------------------------
The BSMestim command has the following parameters and options:
series y: the univariate series whose BSM is to be estimated
series *BSMtrend: name under which to save the trend
series *BSMslope: name under which to save the slope
series *BSMseas: name under which to save the seasonal component
bool restrict_irreg: 0/1 [default=0] restrict irregular variance to zero?
bool restrict_level: 0/1 [default=0] restrict level variance to zero?
bool restrict_slope: 0/1 [default=0] restrict slope variance to zero?
bool restrict_seas: 0/1 [default=0] restrict seasonal variance to zero?
scalar sigmatol: [default=1.E-4] lower bound for zero restriction
If the estimate of any variance goes below sigmatol, the program will restrict
it to zero.
The algorithm starts concentrating the likelihood on sigma^2_e, and runs a
short number of iterations (with a low tolerance) to detect the maximum
variance of the model. In a second step this maximum variance is used to
concentrate the likelihood (using now the default tolerance of gretl).
You can restrict any of the variances of the model to zero by selecting the
corresponding tickmarks in the dialog box.
The q matrix, containing q1, q2 and q3, may be saved assigning a name in the
dialog box.
WARNING: depending on the variance used for concentrating the likelihood, q1,
q2 and q3 may represent different ratios. For example if the model likelihood
is concentrating with respect to sigma^2_xi:
q1 = sigma^2_e/sigma^2_xi,
q2 = sigma^2_eta/sigma^2_xi and
q3 = sigma^2_omega/sigma^2_xi
*/
# private functions
function series kf_filt (series y,
matrix a0,
matrix p0,
matrix bT,
matrix Z,
matrix Sigma_w,
scalar sigma_e,
matrix *at,
matrix *pstar,
series *V,
series *F,
scalar *logLc)
/*
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 i=1..T
# 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
endloop
# 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"
return filtered
end function
function series kf_smooth (matrix pstar,
matrix *at,
matrix bT)
/*
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 i=1..T1
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])
endloop
series ret = at[1,]
#printf "\nSmoothing done\n"
return ret
end function
function scalar BSM (matrix *param,
series y,
matrix *bT,
matrix *at,
matrix *pstar,
series *V,
series *F,
scalar *logLc,
scalar sigmatol[0.0001],
scalar concent[1],
matrix fixed)
if concent > 4
funcerr "concent must be 1, 2, 3, or 4"
endif
genr time
scalar s = $pd
scalar sstart = int(min(time))
scalar send = int(max(time))
scalar T = send-sstart+1
matrix b1 = { 1, 1; 0, 1 }
matrix b2 = zeros(2, s-1)
matrix b3 = -ones(1, s-1) | (I(s-2) ~ zeros(s-2, 1))
matrix bT = (b1 ~ b2) | (b2' ~ b3)
scalar m = cols(bT)
matrix Z = ones(T,1) ~ zeros(T,1) ~ ones(T,1) ~ zeros(T, s-2)
matrix a0 = y[1] * ones(m,1)
matrix p0 = 400000*I(m)
#
matrix Sigma_w = zeros(s+1,s+1)
if concent==1
scalar sigma_e=1
scalar tmp = exp(2*param[1])*fixed[2]
Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
param[1] = (tmp>sigmatol) ? param[1] : -500
scalar tmp = exp(2*param[2])*fixed[3]
Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
param[2] = (tmp>sigmatol) ? param[2] : -500
scalar tmp = exp(2*param[3])*fixed[4]
Sigma_w[3,3] = (tmp>sigmatol) ? tmp : 0
param[3] = (tmp>sigmatol) ? param[3] : -500
elif concent==2
scalar tmp = exp(2*param[1])*fixed[1]
scalar sigma_e=(tmp>sigmatol) ? tmp : 0
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
param[2] = (tmp>sigmatol) ? param[2] : -500
scalar tmp = exp(2*param[3])*fixed[4]
Sigma_w[3,3] = (tmp>sigmatol) ? tmp : 0
param[3] = (tmp>sigmatol) ? param[3] : -500
elif concent==3
scalar tmp = exp(2*param[1])*fixed[1]
scalar sigma_e=(tmp>sigmatol) ? tmp : 0
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
param[2] = (tmp>sigmatol) ? param[2] : -500
scalar tmp = exp(2*param[3])*fixed[4]
Sigma_w[3,3] = (tmp>sigmatol) ? tmp : 0
param[3] = (tmp>sigmatol) ? param[3] : -500
elif concent==4
scalar tmp = exp(2*param[1])*fixed[1]
scalar sigma_e=(tmp>sigmatol) ? tmp : 0
param[1] = (tmp>sigmatol) ? param[1] : -500
Sigma_w[3,3] = 1
scalar tmp = exp(2*param[2])*fixed[2]
Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
param[2] = (tmp>sigmatol) ? param[2] : -500
scalar tmp = exp(2*param[3])*fixed[3]
Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
param[3] = (tmp>sigmatol) ? param[3] : -500
endif
kf_filt(y, a0, p0, bT, Z, Sigma_w, sigma_e, &at, &pstar, &V, &F, &logLc)
return logLc
end function
# public functions
function list BSMestim (series y,
matrix *q[null],
bool irreg[1],
bool level[1],
bool slope[1],
bool seas[1],
scalar sigmatol[0.0001])
set messages off
if (irreg!=1&&irreg!=0) || (level!=1&&level!=0) || (slope!=1&&slope!=0) || (seas!=1&&seas!=0)
funcerr "irreg, level, slope and seas must be 0 or 1"
endif
if $pd<2
funcerr "BSMestim error: your data are not seasonal, you should use the LLTestim function"
endif
matrix fijos = { irreg, level, slope, seas }
matrix theta = { -0.5, -1.5, -2 }
matrix at = { null }
series V = 0
series F = 0
scalar logLc=0
matrix pstar = { null }
matrix bT = { null }
set bfgs_toler 1.E-07
M = BFGSmax(theta, "BSM(&theta, y, &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
scalar q3=exp(2*theta[3])*sstar
matrix qp = { sstar,q1, q2 , q3 }
matrix conc=imaxr(qp)
scalar concent = conc[1]
set bfgs_toler default
matrix theta = { -0.5, -1.5, -2 }
M = BFGSmax(theta, "BSM(&theta, y, &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])
scalar q3=exp(2*theta[3])
string vn=argname(y)
series @vn = y
series @vn_trend = kf_smooth(pstar, &at, bT)
series @vn_slope = at[2,]
series @vn_seas = at[3,]
series @vn_irr = @vn - @vn_trend - @vn_seas
list compo = @vn_trend @vn_slope @vn_seas @vn_irr
printf "\nBasic Structural Model estimation:\n"
printf "-----------------------------------------\n"
printf " loglikelihood\t%#.8g\n", M
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
printf " q3=%8.5f,\t Var(omega)=%8.6E\n", q3, q3*sstar
matrix q = {1, q1, q2, q3}'
elif 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
printf " q3=%8.5f,\t Var(omega)=%8.6E\n", q3, q3*sstar
matrix q = {q1, 1, q2, q3}'
elif 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
printf "q3=%8.5f,\t Var(omega)=%8.6E\n", q3, q3*sstar
matrix q = {q1, q2, 1, q3}'
else
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 "q3=%8.5f,\t Var(xi)=%8.6E\n", q3, q3*sstar
printf "sigma*=\t\t Var(omega)=%8.6E\n", sstar
matrix q = {q1, q2, q3, 1}'
endif
printf "------------------------------------------\n \n"
return compo
end function
##########################
# Original sample script of the package:
open data9-3.gdt
list listavar = BSMestim(reskwh)
print reskwh listavar -o
|
