*
* Monte Carlo exercises from pages 611-613
*
* Because RATS has subscripts based at 1, to get x(0),...,x(100),
* we need to use 101 data points.
*
all 101
*
* (1) Power analysis
*
compute rho=.95
*
compute ndraws=10000
*
dec rect ss
set trend = t
compute countmu=0.0,counttau=0.0
do draws=1,ndraws
*
*    The SET instruction works fine in RATS for generating an AR(1)
*    process. The "loop" here is within the SET instruction, and so
*    is quite a bit more efficient than an actual RATS (or GAUSS)
*    DO loop.
*
   set(first=0.0) y = rho*y{1}+%ran(1.0)
*
*    This keeps computations to a minimum by computing a cross product
*    matrix and then using %sweep to do the regressions, since all the
*    information we need out of the regressions will be included in
*    the results from %sweep.
*
   cmom
   # y y{1} constant trend
   compute ss=%cmom
*
*    Regress on variables 2 (lag of y) and 3 (constant)
*
   compute ss=%sweep(ss,2)
   compute ss=%sweep(ss,3)
*
*    The coefficients of the regression will be in (2,1) and (3,1).
*    (1,1) will be the sum of squared residuals, and (2,2) will be
*    the (x'x)**-1 element for y lagged. The t-stat can be constructed
*    from these pieces of information.
*
   compute tmu=(ss(2,1)-1)/sqrt(ss(1,1)*ss(2,2)/(%nobs-2))
*
*    Now sweep out the trend variable as well
*
   compute ss=%sweep(ss,4)
   compute ttau=(ss(2,1)-1)/sqrt(ss(1,1)*ss(2,2)/(%nobs-3))
   compute countmu=countmu+(tmu<-2.86)
   compute counttau=counttau+(ttau<-3.41)
end do draws
*
disp "Power for t-mu" countmu/ndraws
disp "Power for t-tau" counttau/ndraws
*
* (2) Size distortions
*
compute ndraws=10000
compute theta=-.8
compute countt=0.0,countp=0.0
do draws=1,ndraws
*
*     Because there's only one regression, and not two, the
*     time advantage in avoiding linreg is much lower, so
*     we do the more straightforward coding. The sum of the
*     lag coefficients on the differenced series is obtained
*     using the functions %sum (computes the sum of the entries
*     of a vector or matrix) and %xsubvec (pulls a subvector
*     out of a larger vector, in this case, taking a block of
*     coefficients out of the %beta vector).
*
   set(first=0.0) e = %ran(1.0)
   set(first=0.0) y = y{1}+e+theta*e{1}
   set dy = y-y{1}
   linreg(noprint) y
   # constant y{1} dy{1 to 4}
   compute tmu=(%beta(2)-1)/%stderrs(2)
   compute pmu=%nobs*(%beta(2)-1)/(1-%sum(%xsubvec(%beta,3,6)))
   compute countt=countt+(tmu<-2.86)
   compute countp=countp+(pmu<-14.1)
end do draws
disp "Size for t-mu" countt/ndraws
disp "Size for p-mu" countp/ndraws