*
* Replication file for Baillie and Bollerslev, "The Message in Daily
* Exchange Rates: A Conditional Variance Tale", JBES 1989, vol 7, pp
* 297-305
*
open data bbjbes1989.xls
data(format=xls,org=columns) 1 1245 uk wg ja ca fr it sw dow d81
*
labels fr it ja sw uk wg
# "France" "Italy" "Japan" "Switzerland" "U.K." "West Germany"
*
* Raw data are exchange rates vs the US Dollar
*
dofor s = fr it ja sw uk wg
   log s
end dofor s
*
* Unit root tests
* Table 1
*
report(action=define,title="Table 1. Phillips-Perron Tests for a Unit Root")
dofor s = fr it ja sw uk wg
   report(col=new,atrow=1) %l(s)
   @ppunit(lags=22,det=trend) s
   report(col=current,atrow=2) %cdstat
   @ppunit(lags=22) s
   report(col=current,atrow=3) %cdstat
end dofor
report(action=format,picture="*.###")
report(action=show)
*
* Table 2 estimates
*
report(action=define,title="Daily Normal Distributions")
dofor s = fr it ja sw uk wg
   report(col=new,atrow=1) %l(s)
   clear dx
   set dx = 100.0*(s{0}-s{1})
   linreg(noprint) dx
   # constant
   report(col=curr,atrow=2) %beta(1)
   report(col=curr,atrow=3,special=parens) %stderrs(1)
   report(col=curr,atrow=4) %sigmasq
   report(col=curr,atrow=5) %logl
   set ustd = %resids/sqrt(%seesq)
   corr(print,qstats,number=15,method=yule) ustd
   report(col=curr,atrow=6) %qstat
   set usqr = ustd^2
   corr(print,qstats,number=15,method=yule) usqr
   report(col=curr,atrow=7) %qstat
   stats(noprint) %resids
   report(col=curr,atrow=8) %skewness
   report(col=curr,atrow=9) %kurtosis
end dofor s
report(action=format,picture="*.###",align=decimal)
report(action=show)
*
* Table 3
*
* The DOW series is 2 for Monday,...,6 for Friday. Create dummies for
* each day of the week.
*
dec vect[series] dd(5)
do i=1,5
   set dd(i) = dow==(i+1)
end do i
*
* D81 is a dummy which is 1 for dates before October 1, 1981. The W81
* dummy are the Wednesdays for which that is 1.
*
set w81 = d81.and.dow==4
*
* SKIP is a dummy variable which will be 1 if and only if there is a
* skipped weekday.
*
set(first=1.0) skip = .not.(dow{1}+1==dow.or.(dow{1}==6.and.dow==2))
*
nonlin(parmset=meanparms) b0 b1 b2 b3 b4 b5=-(b1+b2+b3+b4) b8
nonlin(parmset=garchshifts) w0 w1 w2 w3 w4 w5=-(w1+w2+w3+w4) w6 w7
nonlin(parmset=garchparms) alpha1 beta1 rnu
declare series uu h u
frml hmeanf = w0+w1*dd(1)+w2*dd(2)+w3*dd(3)+w4*dd(4)+w5*dd(5)+w6*skip+w7*skip{1}
frml varf = alpha1*uu{1}+beta1*h{1}+hmeanf-(alpha1+beta1)*hmeanf{1}
frml meanf = b0+b1*dd(1)+b2*dd(2)+b3*dd(3)+b4*dd(4)+b5*dd(5)+b8*w81
frml logl = (u=dx-meanf),(uu(t)=u^2),(h(t)=varf(t)),%logtdensity(h,u,1./rnu)
*
* Create two parallel reports, one for the model estimates, one for the summary statistics
*
report(action=define,use=garchestimates)
report(action=define,use=garchsummary,hlabels=||"","France","Italy","Japan","Switzerland","U.K.","West Germany"||)
report(use=garchsummary,atrow=1,atcol=1,fillby=columns) "$Log L$" "$Q(15)$" "$Q^2(15)$" "$m_3$" "$m_4$" $
   "$3(\hat \nu  - 2)(\hat \nu  - 4)^{ - 1}$"
dofor s = fr it ja sw uk wg
   clear dx
   set dx = 100.0*(s{0}-s{1})
   *
   * Get preliminary guess values for the mean parameters estimating just the mean model
   *
   nlls(parmset=meanparms,frml=meanf) dx
   *
   * Use the variance of the estimation to initialize the uu and h series
   *
   compute hinit=%seesq
   *
   set uu = %resids^2
   set h  = hinit
   set u  = %resids
   *
   * Estimation of the full GARCH model will use 1 less data point than
   * the mean model.
   *
   compute gstart=%regstart()+1,gend=%regend()
   *
   * Use base guess values for the variance mean function to avoid
   * problems with negative variances.
   *
   compute w0=hinit,w1=w2=w3=w4=w5=w6=w7=0.0
   *
   * Set initial guess values for the reciprocal degrees of freedom, and
   * the GARCH parameters.
   *
   compute rnu=.10
   compute alpha1=.1,beta1=.8
   *
   maximize(parmset=meanparms+garchshifts+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   report(use=garchestimates,regressors,extra=stderrs)
   report(use=garchsummary,col=new,atrow=1) %funcval
   stats u gstart gend
   set ustd gstart gend = u/sqrt(h)
   corr(print,qstats,number=15,method=yule) ustd gstart gend
   report(use=garchsummary,col=current,atrow=2) %qstat
   set usqr gstart gend = ustd^2
   corr(print,qstats,number=15,method=yule) usqr gstart gend
   report(use=garchsummary,col=current,atrow=3) %qstat
   stats(noprint) ustd
   report(use=garchsummary,col=current,atrow=4) %skewness
   report(use=garchsummary,col=current,atrow=5) %kurtosis
   report(use=garchsummary,col=current,atrow=6) 3.0*(1.0/rnu-2.0)/(1.0/rnu-4.0)
end dofor
report(action=define,use=Table3,title="Table 3 Daily GARCH Models",$
    hlabels=||"","France","Italy","Japan","Switzerland","U.K.","West Germany"||)
report(use=Table3,atcol=1) garchestimates
report(use=Table3,atcol=1) garchsummary
report(use=Table3,action=show)
*
* Table 4
*
report(action=define,title="Likelihood Ratio Tests for Daily Models")
report(atrow=1,atcol=1) "Tests" "France" "Italy" "Japan" "Switzerland" "U.K" "West Germany"
dofor s = fr it ja sw uk wg
   clear dx
   set dx = 100.0*(s{0}-s{1})
   *
   * Get preliminary guess values for the mean parameters estimating just the mean model
   *
   nlls(parmset=meanparms,frml=meanf) dx
   *
   * Use the variance of the estimation to initialize the uu and h series
   *
   compute hinit=%seesq
   *
   set uu = %resids^2
   set h  = hinit
   set u  = %resids
   *
   * Estimation of the full GARCH model will use 1 less data point than
   * the mean model.
   *
   compute gstart=%regstart()+1,gend=%regend()
   compute w0=hinit,w1=w2=w3=w4=w5=w6=w7=0.0
   *
   * Set initial guess values for the reciprocal degrees of freedom, and
   * the GARCH parameters.
   *
   compute rnu=.10
   compute alpha1=.1,beta1=.8
   *
   * Alternative PARMSETS for doing the various restrictions
   *
   nonlin(parmset=meanrestricted) b0 b1=0 b2=0 b3=0 b4=0 b5=0 b8=0
   nonlin(parmset=varrestricted) w0 w1=0 w2=0 w3=0 w4=0 w5=0 w6=0 w7=0
   nonlin(parmset=varwoholidays) w0 w1 w2 w3 w4 w5=-(w1+w2+w3+w4) w6=0 w7=0
   nonlin(parmset=nogarch) alpha1=0 beta1=0 rnu
   *
   * Base model with all parameters unrestricted
   *
   maximize(parmset=meanparms+garchshifts+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute baselogl=%logl,baseregs=%nreg
   *
   * Without any shifts in the mean
   *
   maximize(parmset=meanrestricted+garchshifts+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute teststat=2*(baselogl-%logl),signif=%chisqr(teststat,baseregs-%nreg)
   report(atrow=2,col=new) teststat
   report(atrow=3,col=current,special=parens) signif
   maximize(parmset=meanparms+varwoholidays+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute teststat=2*(baselogl-%logl),signif=%chisqr(teststat,baseregs-%nreg)
   report(atrow=4,col=current) teststat
   report(atrow=5,col=current,special=parens) signif
   maximize(parmset=meanparms+varrestricted+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute teststat=2*(baselogl-%logl),signif=%chisqr(teststat,baseregs-%nreg)
   report(atrow=6,col=current) teststat
   report(atrow=7,col=current,special=parens) signif
   maximize(parmset=meanrestricted+varrestricted+garchparms,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute teststat=2*(baselogl-%logl),signif=%chisqr(teststat,baseregs-%nreg)
   report(atrow=8,col=current) teststat
   report(atrow=9,col=current,special=parens) signif
   maximize(parmset=meanrestricted+varrestricted+nogarch,pmethod=simplex,piter=10,method=bhhh) logl gstart gend
   compute teststat=2*(baselogl-%logl),signif=%chisqr(teststat,baseregs-%nreg)
   report(atrow=10,col=current) teststat
   report(atrow=11,col=current,special=parens) signif
end dofor
report(atcol=1,atrow=2,fillby=columns) "Mean Shifts" "" "Holiday Shifts" "" "Variance Shifts" "" $
    "Mean+" " Variance Shifts" "Mean+Variance+" " GARCH coefficients"
report(action=format,align=decimal,picture="*.###")
report(action=show)
*
* Table 5
* Wednesdays only (or Thursday if Wednesday is missing)
*
set pullweds = dow==4.or.dow==5.and.dow{1}==3
dofor s = fr it ja sw uk wg
   sample(smpl=pullweds) s / wx
   clear dx
   set dx = 100.0*(wx-wx{1})
   garch(p=1,q=1) / dx
end dofor s
*
* Table 6
* Every two week data
*
dofor s = fr it ja sw uk wg
   sample(smpl=pullweds) s / wx
   sample(smpl=%clock(t,2)==1) wx / w2x
   clear dx
   set dx = 100.0*(w2x-w2x{1})
   garch(q=1) / dx
end dofor s
*
* Table 7
*
dofor s = fr it ja sw uk wg
   sample(smpl=pullweds) s / wx
   sample(smpl=%clock(t,4)==1) wx / w4x
   clear dx
   set dx = 100.0*(w4x-w4x{1})
   garch(q=0,p=0) / dx
end dofor s