State space applications
The following program estimates an ARIMA(2,1,2) model with NLS and
with full maximum likelihood (using the state-space object). To find the
state-space representation of an unrestricted ARMA model, use the autospec
menu interactively. This program calls a subroutine to get starting values
for the unrestricted ARMA parameters using the fast Hannan-Rissanen algorithm.
'estimation of ARIMA(2,1,2) model
'version 4
'9/21/2000 h
include sub_HannanRissanen.prg
workfile arima22 q 1947:1 1986:3
series y
y.fill 1056.5, 1063.2, 1067.1, 1080.0, 1086.8, 1106.1, 1116.3, 1125.5, 1112.4, 1105.9, 1114.3, 1103.3, 1148.2, 1181.0, 1225.3, 1260.2, 1286.6, 1320.4, 1349.8, 1356.0, 1369.2, 1365.9, 1378.2, 1406.8, 1431.4, 1444.9, 1438.2, 1426.6, 1406.8, 1401.2, 1418.0, 1438.8, 1469.6, 1485.7, 1505.5, 1518.7, 1515.7, 1522.6, 1523.7, 1540.6, 1553.3, 1552.4, 1561.5, 1537.3, 1506.1, 1514.2, 1550.0, 1586.7, 1606.4, 1637.0, 1629.5, 1643.4, 1671.6, 1666.8, 1668.4, 1654.1, 1671.3, 1692.1, 1716.3, 1754.9, 1777.9, 1796.4, 1813.1, 1810.1, 1834.6, 1860.0, 1892.5, 1906.1, 1948.7, 1965.4, 1985.2, 1993.7, 2036.9, 2066.4, 2099.3, 2147.6, 2190.1, 2195.8, 2218.3, 2229.2, 2241.8, 2255.2, 2287.7, 2300.6, 2327.3, 2366.9, 2385.3, 2383.0, 2416.5, 2419.8, 2433.2, 2423.5, 2408.6, 2406.5, 2435.8, 2413.8, 2478.6, 2478.4, 2491.1, 2491.0, 2545.6, 2595.1, 2622.1, 2671.3, 2734.0, 2741.0, 2738.3, 2762.8, 2747.4, 2755.2, 2719.3, 2695.4, 2642.7, 2669.6, 2714.9, 2752.7, 2804.4, 2816.9, 2828.6, 2856.8, 2896.0, 2942.7, 3001.8, 2994.1, 3020.5, 3115.9, 3142.6, 3181.6, 3181.7, 3178.7, 3207.4, 3201.3, 3233.4, 3157.0, 3159.1, 3199.2, 3261.1, 3250.2, 3264.6, 3219.0, 3170.4, 3179.9, 3154.5, 3159.3, 3186.6, 3258.3, 3306.4, 3365.1, 3444.7, 3487.1, 3507.4, 3520.4, 3547.0, 3567.6, 3603.8, 3622.3, 3655.9, 3661.4, 3683.3
'take log 1st difference
series dlogy = dlog(y)
'NLS/conditional MLE
equation eq1.ls(showopts,c=1e-5,m=1000) dlogy c ar(1) ar(2) ma(1) ma(2)
show eq1.output
'setup ARMA(2,2) in sspace
sspace ss1
ss1.append @signal dlogy = c(1) + sv1 + c(2)*sv2 + c(3)*sv3
ss1.append @state sv1 = c(5)*sv1(-1) + c(6)*sv2(-1) + [var = exp(c(4))]
ss1.append @state sv2 = sv1(-1)
ss1.append @state sv3 = sv2(-1)
'these starting values (same as NLS) take 64 iterations to converge
'c(1)=0.00794 'constant
'c(2)=0.0025 'MA(1)
'c(3)=0.0025 'MA(2)
'c(4)=log(@var(dlogy)) 'sigma
'c(5)=0.0025 'AR(1)
'c(6)=0.0025 'AR(2)
'starting values from Hannan-Rissanen
vector(6) p0
call sub_HannanRissanen(dlogy, 2, 2, 8, p0)
c(1)=p0(1) 'constant
c(2)=p0(4) 'MA(1)
c(3)=p0(5) 'MA(2)
c(4)=log(p0(6)) 'sigma
c(5)=p0(2) 'AR(1)
c(6)=p0(3) 'AR(2)
'estimate
ss1.ml(showopts,m=1000,c=1e-5)
freeze(out1) ss1.output
show out1
'different starting values which fail to improve
c(1)=0 'constant
c(2)=-1.328 'MA(1)
c(3)=0.392 'MA(2)
c(4)=log(p0(6)) 'sigma
c(5)=1.68 'AR(1)
c(6)=-0.749 'AR(2)
ss1.ml(showopts,m=1000,c=1e-5)
freeze(out2) ss1.output
show out2
The following program estimates the multiplicative seasonal ARMA model, known
as the "airline" model, using NLS and full MLE. This example shows how to specify
a multiplicative seasonal MA model in state-space form.
'estimation of airline model
'version 4
'9/21/2000 h
workfile airline m 1949:01 1960:12
series x
x.fill 112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118, 115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140, 145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166, 171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194, 196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201, 204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229, 242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278, 284, 277, 317, 313, 318, 374, 413, 405, 355, 306, 271, 306, 315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336, 340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337, 360, 342, 406, 396, 420, 472, 548, 559, 463, 407, 362, 405, 417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432
'airline model transformation
series y = dlog(x,1,12)
'NLS/conditional MLE
equation eq1.ls(showopts,c=1e-5,m=1000) y ma(1) sma(12)
show eq1.output
'setup ARMA(0,1)(0,1) in sspace
sspace ss1
ss1.append @signal y = sv1 + c(1)*sv2 + c(2)*sv13 + c(1)*c(2)*sv14
ss1.append @state sv1 = [var = exp(c(3))]
ss1.append @state sv2 = sv1(-1)
ss1.append @state sv3 = sv2(-1)
ss1.append @state sv4 = sv3(-1)
ss1.append @state sv5 = sv4(-1)
ss1.append @state sv6 = sv5(-1)
ss1.append @state sv7 = sv6(-1)
ss1.append @state sv8 = sv7(-1)
ss1.append @state sv9 = sv8(-1)
ss1.append @state sv10 = sv9(-1)
ss1.append @state sv11 = sv10(-1)
ss1.append @state sv12 = sv11(-1)
ss1.append @state sv13 = sv12(-1)
ss1.append @state sv14 = sv13(-1)
'these starting values (same as NLS) do not converge after 1000 iterations
'c(1)=0.0025 'MA(1)
'c(2)=0.01 'MA(12)
'c(3)=log(@var(y)) 'sigma
'use NLS estimates as starting values
c(3)=log(@var(y)) 'sigma
'estimate with BHHH (Marquardt does not converge after 1000 iterations)
ss1.ml(showopts,m=1000,c=1e-5,b)
show ss1.output
The following program illustrates the use of the state-space object as
a filtering/smoothing routine without estimating any parameters.
'cubic spline smooothing using sspace
'last checked 10/27/2000 h
'change path to program path
%path = @runpath
cd "{%path}"
' load workfile
load nber
series y = m03054
smpl @first 42:06
'set smoothing parameter (signal-to-noise ratio)
!q = 0.005
'setup sspace
sspace ss1
ss1.append @signal y = sv1 + [var=1.0]
ss1.append @state sv1 = sv1(-1) + sv2(-1) + [ename=u1]
ss1.append @state sv2 = sv2(-1) + [ename=u2]
'c(1)*c(1) is the smoothing parameter (signal-to-noise ratio)
'ss1.append @evar cov(u1,u1)=c(1)*c(1)/3
'ss1.append @evar cov(u1,u2)=0.5*c(1)*c(1)
'ss1.append @evar cov(u2,u2)=c(1)*c(1)
ss1.append @evar cov(u1,u1)=!q/3
ss1.append @evar cov(u1,u2)=0.5*!q
ss1.append @evar cov(u2,u2)=!q
'parse (no parameters to estimate)
ss1.ml(d,m=100,c=1e-9)
'get smoothed estimate
ss1.makesignals(t=smooth) ys1
'get HP filtered series
!lambda = 14400
y.hpf(!lambda) ys2
'plot
graph graph1.line y ys1 ys2
graph1.setelem(1) legend(Actual)
graph1.setelem(2) legend(Cubic spline smoothing (q = !q))
graph1.setelem(3) legend(Hodrick-Prescott filter (lambda = !lambda))
graph1.legend position(.5,0) -inbox
graph1.options size(6,2)
show graph1
The following program illustrates how to estimate a basic stochastic
volatility model by quasi maximum likelihood (QML). The program compares
and plots the estimated log volatility series to that from a GARCH(1,1)
model. Note that the standard errors reported in the state space estimates
are not correct (we currently do not have an option to compute QML standard
errors in statespace).
'basic stochastic volatility model
'QML estimation by state space
'12/14/2000 h.
create u 1 1000
'get data
%datafile = @runpath + "svpdx.dat"
read(skiprow=1) "{%datafile}" rt
'estimate GARCH(1,1) model
equation eq1.arch(1,1,m=500,c=1e-5,showopts) rt c
show eq1.output
'and get conditional variance series
smpl 946 1000
eq1.forecast rhat rhat_se gt
smpl @all
'transformation for QML estimation
series y = log(rt*rt)
'variance of log chi-square distribution
!pi = @acos(-1)
scalar s2 = 0.5*!pi*!pi
'specify quasi-likelihood state-space
sspace sv
sv.append @signal y = -1.27 + ht + [var=s2]
sv.append @state ht = c(1) + c(2)*ht(-1) + [var=exp(c(3))]
'set starting values
c(1) = 0.01
c(2) = 0.9
c(3) = 0.1
'estimate
show sv.ml(showopts,m=500,c=1e-7)
'and retrieve state series
sv.makestate(t=pred) hf
sv.makestate(t=predse) hf_se
sv.makestate(t=filt) ht
sv.makestate(t=filtse) ht_se
sv.makestate(t=smooth) hs
sv.makestate(t=smoothse) hs_se
'plot to compare log variance from GARCH and stochastic volatility
graph graph1.line log(gt) hf hs
graph1.options size(8,2)
graph1.addtext(0.1,-0.3) Log volatility
graph1.setelem(1) legend(GARCH(1,1))
graph1.setelem(2) legend(One-step ahead)
graph1.setelem(3) legend(Smoothed)
graph1.legend columns(1) position(0,0)
show graph1