sspace

# State space applications

- ARIMA(2,1,2) (arima22.prg)
- Seasonal ARIMA(0,1,1)(0,1,1) airline model (airline.prg)
- Cubic spline smoothing (spline.prg)
- QML estimation of stochastic volatility (svol.prg)

## ARIMA(2,1,2)

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

## Seasonal ARIMA(0,1,1)(0,1,1) (airline model)

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

## Cubic spline smoothing

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

## QML estimation of stochastic volatility

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