options crt;  ? 2e10  rev. 9/97
in return;
freq m;
smpl 78:1 87:12;
rpmark = market - rkfree;
rpmk2 = rpmark*rpmark;
dot mobil texaco ibm dec datgen coned psnh weyer boise
    ?motor  ? skip Motorola for now, due to its strange sample
    tandy panam delta contil citcrp gerber genmil;
  rp. = . - rkfree;
  regopt(pvprint,lmlag=1,dwpval=exact) whiteht,lmar;  ? turn on tests
  olsq rp. c rpmark;
  unmake @coef alpha beta;  ? save coefficients for MA(1) starting values
  rename @whiteht whiteht;
  rename @lmar1 lmar1;
  rename @dw dwstat;
  rename %dw dwpval;
  rename @res e;      ? save residuals and DW-stat for later
  ? turn off these extra tests for regressions below, for less confusion
  regopt(nocalc,lmlag=0,dwpval=approx) whiteht,lmar;
  ? calculate regression diagnostics by hand for comparison
  e2 = e*e;
  olsq(silent) e2 c rpmark rpmk2;
  set hetstat = @rsq*@nob;
  print hetstat whiteht;  ? print both to confirm they are equal
  cdf(chisq,df=2,print) hetstat pvhet;
  if pvhet < .05; then;
   do;
    title 'Reject homoscedasticity';
    wt = 1/@fit;
    olsq(weight=wt) rp. c rpmark;
   endd;
  else;
    title 'Cannot reject homoscedasticity';  ? (a)
?
  print dwstat,dwpval;
  ?if dwstat < 1.69; then;  ? 5%, t=100, k'=1, du critical value
  if dwpval < .05; then;  ? using exact P-value
   do;
    ?title 'Reject rho=0 or inconclusive';  ? for inconclusive version
    title 'Reject rho=0';  ? for exact version
    ar1 rp. c rpmark;  ? improved Cochran-Orcutt method (ML)
   endd;
  else;
    title 'Cannot reject rho=0';   ? (b)
?
  smpl 78:2 87:12;
   olsq(silent) e c rpmark e(-1);  ? LM test for MA(1) or AR(1)
   set mastat = @t(3);
   set mastatpf = 2*@fst;  ? @LMAR uses p*F version of LM test
   print mastat mastatpf lmar1;  ? print both to confirm they are equal
   set dft = @nob-@ncid;
   cdf(t,df=dft,print) mastat pmastat;
   if pmastat >= .05; then;
    title 'Cannot reject MA(1)=0 - LM';   ? (c)
   else; do;
    title 'Reject MA(1)=0 - LM';
?
?    Note:  regression estimation with MA(1) error term involves a lot
?    of programming at present (5/95).  TSP 4.3 commands are used.
?----------------------------------------------------
? ARMAX(0,1) estimation
? following Andrew Harvey, The Econometric Analysis of Time Series 1990 p.273
?
  param theta alpha beta;
  smpl 78:1 87:12;
  eps = 0; deda = 0; dedb = 0; dedt = 0;
  smpl 78:2 87:12;
  eps = e;
?
? obtain initial value for theta (MA(1) parameter)
?
  bjest(nma=1,nocumplo,noplot,silent) eps;
  set theta = -@coef(1);  ? sign of coef is reversed in BJEST
  list bnames alpha beta theta;
  mmake b bnames;
?
? iteration
?
? residual
  frml eqe eps = rp. - alpha - beta*rpmark - theta*eps(-1);
? derivatives
  frml d1 deda = -1       - theta*deda(-1);
  frml d2 dedb = -rpmark  - theta*dedb(-1);
  frml d3 dedt = -eps(-1) - theta*dedt(-1);
?
  genr(silent) eqe;
  mat ssrold = eps'eps;
?
  mat novar = nrow(b);
  ? Turn off diagnostics for this regression, since we only
  ? want the coefficients.
  supres @logl @nob @coef;
  regopt(nocalc,resetord=0) dw lmhet jb fst;
  const maxit,20 maxsqz,10 tol,.001;
do iter=1,maxit;
  genr(silent) d1; genr(silent) d2; genr(silent) d3;
  ? Gauss-Newton iteration:  regress eps on derivatives
  ols(silent) eps deda dedb dedt;
  rename @coef delt;
  ? test for convergence
  mat worst = max(abs(delt)/(abs(b) + 10*tol));
  if worst<tol; then; goto 20;
  do isqz=1,maxsqz;
    mat testb = b - delt/(4**(isqz-1));
    unmake testb bnames;
    genr(silent) eqe;
    mat ssrnew = eps'eps;
    ?write ssrold ssrnew isqz;
    if ssrnew<ssrold; then; goto 15;
  enddo;
  write(form="(' Failure to improve after ',F4.0,' squeezes.')") isqz;
  write b delt;
  goto 18;
  15 ?write iter,ssrnew;
  set ssrold = ssrnew;
  copy testb b;
enddo;
18 write(form="(' Failure to converge after ',F4.0,' iterations.')") iter;
goto 30;
20 write(form="(' Convergence after ',F4.0,' iterations.')") iter;
30 unmake b bnames;
?write ssrold bnames;
nosupres;
tstats(names=bnames) b @vcov;  ? (c)  A t-test may be used on theta
  if abs(@t(3)) > 1.96; then;  ? 5% critical value for asymptotic normal
    title 'Reject MA(1)=0 - Wald';
   else;
    title 'Cannot reject MA(1)=0 - Wald';   ? (c)
 enddo;
?----------------------------------------------------
?    End of regression with MA(1) estimation
?
  smpl 78:1 87:12;
?
? Normality tests
?
? 1. Done via regression diagnostics
  regopt(pvprint) jb swilk;  ? Jarque-Bera and Shapiro-Wilk
  olsq e c;   ? (d)
  regopt(nocalc) jb swilk;  ? turn back off
?
? 2. Done by hand -- Jarque-Bera, aka chi-squared
  msd e;
  ? The first formula below for Jarque-Bera is based on the discussion
  ? in the TSP Reference Manual under MSD.  It does not exactly
  ? match @JB from OLSQ, because @SKEW and @KURT contain some
  ? small-sample adjustments.
  set jbstat = @nob*( @skew**2/6 + @kurt**2/24 );
  ? The next formulas below for Jarque-Bera match @JB exactly,
  ? because they take account of the small-sample adjustments
  ? to @SKEW and @KURT to extract m2, m3, and m4.
  set m2 =  @var*(@nob-1)/@nob;
  set m3 =  @skew*(@stddev**3)*(@nob-1)*(@nob-2)/@nob**2;
  set m4 = [@kurt*(@stddev**4)*(@nob-1)*(@nob-2)*(@nob-3)/@nob**2
            + 3*(@nob-1)*m2**2]/(@nob+1);
  set b23 = m4/m2**2 - 3;  ? @kurt without small-sample adjustment
  set jbsame = @nob*( (m3**2/m2**3)/6 + b23**2/24 );
  print jbstat jbsame @jb;
  if jbstat > 4.34; then;  ? J-B 5% critical value for N=125
    title 'Reject normality at 95% level';
  else;
    title 'Cannot reject normality at 95% level';  ? (d)
enddot;