*ARCH_WPI.PRG
cal 1960 1 4             ;* set calendar to begin 1960:Q1
all 0 1992:2             ;* data runs until 1992:Q2
open data wpi.wk1      ;* data is assumed to be on drive a:\
data(format=wks,org=obs)
set dlwpi = log(wpi) - log(wpi{1}) ;* transform data

boxjenk(constant,ar=1,ma=||1,4||) dlwpi / resids

cor(partial=pacf,qstats,number=36,span=12,dfc=4) resids

set ressq = resids*resids     ;* Form the squared residuals
cor(partial=pacf,qstats,number=24,span=4,dfc=4) ressq

linreg ressq             ;* Perform the Lagrange multiplier test for ARCH(4) errors
# constant ressq{1 to 4}

compute trsq = %nobs*%rsquared     ;* Compute T*R2 and obtain the cumulative density
cdf chisqr trsq 4                            ;* of trsq as chi-square with 4 degrees of freedom.

nonlin b0 b1 a0 a1
frml e = dlwpi - b0 - b1*dlwpi{1}
frml var = a0 + a1*e(t-4)**2
frml L = -0.5*(log(var)+e(t)**2/var)
boxjenk(noprint,constant,ar=1) dlwpi
compute b0=%beta(1), b1=%beta(2)
compute a0=%seesq, a1=.8
maximize(iterations=75) L 10 *


set u = 0.0
nonlin b0 b1 b2 a0 a1
frml e = dlwpi - b0 - b1*dlwpi{1} - b2*u{1}
frml var = a0 + a1*e(t-4)**2
frml L = (u = e), -0.5*(log(var)+e(t)**2/var)
boxjenk(noprint,constant,ar=1,ma=1) dlwpi
compute b0=%beta(1), b1=%beta(2), b2 = %beta(3)
compute a0=%seesq, a1=.8
maximize(iterations=75) L 10 *

set u = 0.0
nonlin b0 b1 b2 b3 a0 a1
frml e = dlwpi - b0 - b1*dlwpi{1} - b2*u{1} - b3*u{4}
frml var = a0 + a1*e(t-4)**2
frml L = (u = e), -0.5*(log(var)+e(t)**2/var)
boxjenk(noprint,constant,ar=1,ma=||1,4||) dlwpi
compute b0=%beta(1), b1=%beta(2), b2 = %beta(3) , b3 = %beta(4)
compute a0=%seesq, a1=.8
maximize(iterations=75) L 10 *

set w = 0.0
set u = 0.0
nonlin b0 b1 b2 b3 a0 a1 a2
frml e = dlwpi - b0 - b1*dlwpi{1} - b2*u{1} - b3*u{4}
frml var = a0 + a1*e{1}**2 + a2*w{1}
frml L = (u = e), (w = var), -.5*(log(var)+e(t)**2/var)
boxjenk(noprint,constant,ar=1,ma=||1,4||) dlwpi
compute b0=%beta(1), b1=%beta(2), b2 = %beta(3) , b3 = %beta(4)
compute a0=%seesq, a1=.3 , a2 = .5
maximize(iterations=75) L 10 *


set resids = 0.0
set resids 62:2 92:2 = dlwpi - %beta(1) - %beta(2)*dlwpi{1}  $
   - %beta(3)*resids{1} - %beta(4)*resids{4}
cor(qstats,dfc=4,number=32,span=8) resids 62:2 92:2

source(noecho) bjident.src
@bjident resids 62:2 92:2
set convar = 0.0
set convar 62:2 92:2 = %beta(5) + %beta(6)*resids{1}**2 $
     + %beta(7)*convar{1}
set upper = dlwpi + 2*convar**0.5
set lower = dlwpi - 2*convar**0.5
graph(header='Two Conditional Standard Deviations Around DLWPI',  $
     key=upleft,patterns) 3
# upper  62:2 92:2  1
# lower  62:2 92:2  1
# dlwpi