*
* RATS Program for Chapter 6, empirical exercise 1
* pp 438-440. Analyzes the yen data.
*
cal(w) 1975:1:3
all 778
*
* read in data
*
open data yen.xls
data(org=columns,format=xls) / s f s30

*  The dataset has the following variables in weekly data:
*  Column 1: spot exchange rate, foreign currency price of $
*  Column 2: 30-day forward rate
*  Column 3: spot exchange rate 30 days hence
* From the first week of 1975 to the last week of 1989
*
* Make required transformations
*
set s30_f  = 1200*(log(s30)-log(f))
set s30_s  = 1200*(log(s30)-log(s))
set f_s    = 1200*(log(f)-log(s))
set logs30 = log(s30)
set logf   = log(f)
*
* Graph the forecast error
*
graph(footer="Figure 6.1: Forecast Error, Yen/Dollar",vlabel="s3-f (annualized percentage change)")
# s30_f
*
* Graph the spot rate
*
graph(footer="Figure 6.3: Yen/Dollar Spot Rate, Jan 1975-Dec 1989")
# s
*
* X-Y graph of s30 vs forward
*
scatter(footer="Figure 6.4: Plot of s30 against f, Yen/Dollar",vlabel="s30(t), log spot rate 30 days hence",$
  hlabel="f(t), log forward rate")
# logf logs30

*
* part (b): test for serial correlation
*
* simple statistics
*
table / s30_f
*
* correlogram. The procedure regcorrs is designed for computing and graphing
* a correlogram for diagnostic purposes. The NOCRIT option suppresses the output
* of regression information criteria, since we're not estimating a regression here.
*
@RegCorrs(number=40,nocrit,footer="Figure 6.2 Correlogram of s30-f, Yen/Dollar") s30_f
*
* However, if you want to compute them yourself, use the correlate instruction
*
correlate(qstats,number=40) s30_f
*
* part (c): is the spot rate a random walk?
*
set ds = 1200*(log(s)-log(s{1}))
*
* simple statistics
*
table / ds
*
* correlogram
*
@RegCorrs(number=40,nocrit) ds
*
* part (d): unconditional test. This can be done as a regression on a constant
* with robust standard errors with LAGS=4.
*
linreg(lags=4,lwindow=flat) s30_f
# constant
test(zeros)
# 1
*
* part (e): regression test with truncated kernel
*
linreg(lags=4,lwindow=flat) s30_s
# constant f_s
*
* test H0: b0=0 and b1=1
*
test(all,vector=||0.0,1.0||)
*
* part (f): regression test with Bartlett kernel with lag = 12
*
linreg(lags=12,lwindow=bartlett) s30_s
# constant f_s
*
* test H0: b0=0 and b1=1
*
test(all,vector=||0.0,1.0||)
*
* VARHAC
*
compute pmax=fix(%nobs**.333)
linreg s30_s / resids
# constant f_s
*
set g1 = resids
set g2 = resids*f_s
*
* This makes extensive use of the %SWEEP function, which, when applied to a cross
* product matrix, "sweeps" out the effect of the pivot variable. What remains, at
* the diagonal elements, are the sums of squared residuals from the regression on
* all of the pivots swept out.
*
* The cross product matrix is formed with the two "dependent variables" resids x
* constant and resids x f_s, are in the first two rows, and their lags in the
* remainder.
*
cmom
# g1 g2 g1{1 to pmax} g2{1 to pmax}
*
compute k=2
dec vect[integer] biclags(k)
dec vect bicbest(k)
*
dec rect sxx
compute biclagsall=0
*
* Loop over the lags.
*
compute sxx=%cmom
do p=1,pmax
*
*     Sweep out the current lag for each of the regressors
*
   do j=1,k
      compute jl=k+(j-1)*pmax+p
      compute sxx=%sweep(sxx,jl)
   end do j
*
*    For each regressor, compute its VAR(p) BIC. If this is better than the best
*    value so far for this regressor, save the lag and value. Also, set biclagsall
*    to this value of the lag. biclagsall will, at the end, be equal to the longest
*    lag that was selected for any variable.
*
   do i=1,k
      compute bic=log(sxx(i,i)/%nobs)+p*k*log(%nobs)/%nobs
      disp i p bic
      if p==1.or.bic<bicbest(i) {
         compute bicbest(i)=bic,biclags(i)=p
         compute biclagsall=p
      }
   end do i
end do p
*
* This does the same number of lags for each regressor
*
compute sxx=%cmom
do p=1,biclagsall
   do j=1,k
      compute jl=k+(j-1)*pmax+p
      compute sxx=%sweep(sxx,jl)
   end do j
end do p
*
compute [rect] lagsums=%identity(k)
do p=1,biclagsall
   do i=1,k
      do j=1,k
         compute lagsums(i,j)=lagsums(i,j)-sxx(k+(j-1)*pmax+p,i)
      end do j
   end do i
end do p
*
* Recolor
* The sum of ee' is in the top kxk corner of sxx. We don't divide by the number
* of observations as all the "N" factors will wash out from the calculation being
* done.
*
compute sigma=%xsubmat(sxx,1,k,1,k)
compute shat=%mqform(sigma,tr(inv(lagsums)))
*
linreg(create,lastreg,form=chisquared,covmat=%xx*shat*%xx,title="OLS with VARHAC Std Errors")