@ This GAUSS file reads in data, sets options, and calls numerical
   optimization routine for numerical estimation of SWARCH model @

        output file=junk.out reset;
        format /m1 /ros 16,8;
/* ======================================================================= */
@  Alter this section so as to read in your data and set parameters @

           @ Crisp data @
           capt = 1331;        @ capt is the sample size @
           load y[capt,1] = crisp.w;

@--------- Adjust any of the following to control specification desired --- @
           nk = 5;           @ nk is the first observation to be used in
                               estimation @
           ns = 3;           @ ns is the number of primitive states @
           pphi = 1;         @ pphi is the number of lags in autoregression
                                for y; @
           izz = 1;          @ izz = 1 means pij parameterized as
                                       th(ij)^2/sum(th(ij)^2)
                               izz = 2 means pij parameterized as pij @
           ipm =3;          @ ipm specifies order in which transition probs
                               are parameterized
                                  ipm = 1 implies p11 and p22 estimated
                                  ipm = 2 implies pij for i=1,..,n j=1,..,n-1
                                  ipm = 3 user input code @
           irs = 2;          @ irs denotes the number of probabilities that
                               are restricted @
           karch = 2;        @ karch is the order of the ARCH process @
           larch = 1;        @ larch = 0 means no leverage effect,
                               larch = 1 means leverage effect @
           barch = 2;        @ barch = 0 means estimate initial variance by MLE
                               barch = 2 means don't need initial variance @
           tarch = 1;        @ tarch = 0 means normal
                               tarch = 1 means t distribution @

@------ Input initial values for parameters ------------ @
/*   The order in which variables are represented is as follows
              constant term in regression
              autoregressive terms in regression
              initial variance parameter
              constant term in ARCH equation
              ARCH params in state 1
              next elements:   when izz =0 these are the transition probs
                               when izz =1 these are params v(i,j) such that
                                    p(i,j) = v(i,j)^2 / sum j v(i,j)^2
                               elements are ordered as p11, p22 when ns =2
                                    ordered as p(1,1),p(2,1),...,p(ns,1),
                                    p(1,2),...,p(ns-1,ns) when ns > 2
              next (ns - 1) elements: factor squared residuals are divided by
                                      to get non-switching ARCH
              leverage parameter
              degrees of freedom parameter for t distribution */

      proc startval; @procedure to set starting values @
      local th;
      let th[13,1] =
     0.35080500       0.25003200      -0.56817500
     0.027988000      -0.11706300
     11.4016   0.051132  10.765493  0.1208310
       4.3512660        13.146594       0.42415600       -5.2199370 ;

/*         let th[11,1] =
           0.34775  0.25589  0.55992  0.06036  0.16047  0.051751
           0.99195  0.99907  5.7857  0.44348  3.7547 ; */
      retp(th);
      endp;

/*==========================================================================*/
  @ In general no parts of this section should be changed @

ps = karch;  @ ps is the number of lagged states that matter @
n = ns^(ps+1);     @ n is the dimension of the state vector @
kc = 1;            @ kc = 2 to print out ergodic probs and likelihood value
                     with each call to likelihood @
ks = 1;            @ ks = 2 if smoothed probs are to be calculated @
captst = capt - nk +1; @ captst is the effective sample size @
resid = zeros(capt,1);  @ regression residuals (filled in by procedures)@
skif = zeros(captst,n); @ skif is the matrix of filtered probs @
skis = zeros(captst,n); @ skis is the matrix of smoothed probs @
varfor = zeros(capt,1); @varfor is the vector of forecast variances @
outprob = zeros(capt,1);  @outprob is the vector of outlier probabilities @
id = eye(ns);           @ used in certain calculations below @

proc pattern1; @ This proc returns a (ps+1)*ns x n matrix.  The ith
                column contains a one in row j if st = j, contains a
                one in row ns+j if st-1 = j, and so on @
     local i1,ix,iq,na;
     na = n/ns;
     ix = eye(ns).*.ones(1,na);
     i1 = 1;
     do until i1 > ps;
       na = na/ns;
       iq = ones(1,ns^i1).*.(eye(ns).*.ones(1,na));
       ix = iq|ix;
     i1 = i1+1;
     endo;
retp(ix);
endp;

/* ======================================================================= */
proc matpm(xth);  @This proc defines the user's conventions for reading
                 elements of Markov transition probabilities from
                 parameter vector @
   local pm,ixth;
   ixth = rows(xth);
   pm = zeros(ns,ns);
     if ipm == 1;  @ for ns =2 this option has parameters as p11 and p22 @
          if izz == 1;
             pm[1,1] = xth[1,1]^2/(1 +xth[1,1]^2);
             pm[2,2] = xth[2,1]^2/(1 + xth[2,1]^2);
          else;
             pm[1,1] = xth[1,1];
             pm[2,2] = xth[2,1];
          endif;
          pm[2,1] = 1 - pm[1,1];
          pm[1,2] = 1 - pm[2,2];
     elseif ipm == 2;  @ general case has parameters pij for i = 1,...,n and
                          j = 1,...,n-1 @
        pm[1:ns-1,.] = reshape(xth[1:ixth,1],ns-1,ns);
        if izz == 1;
           pm[ns,.] = ones(1,ns);
           pm = pm^2;
           pm = pm./(sumc(pm)');
        else;
           pm[ns,.] = (1 - sumc(pm))';
        endif;
     elseif ipm == 3;  @ This section can be rewritten by user to impose zeros
                          and ones where desired @
        if izz == 1;
             pm[1,1] = xth[1,1]^2/(1 + xth[1,1]^2);
             pm[1,3] = xth[2,1]^2/(1 + xth[2,1]^2 + xth[4,1]^2);
             pm[2,1] = 1/(1 + xth[1,1]^2);
             pm[2,2] = xth[3,1]^2/(1 + xth[3,1]^2);
             pm[2,3] = xth[4,1]^2/(1 + xth[2,1]^2 + xth[4,1]^2);
             pm[3,2] = 1/(1 + xth[3,1]^2);
             pm[3,3] = 1/(1 + xth[2,1]^2 + xth[4,1]^2);

        elseif izz == 2;
             pm[1,1] = xth[1,1];
             pm[1,3] = xth[2,1];
             pm[2,1] = 1 - xth[1,1];
             pm[2,2] = xth[3,1];
             pm[2,3] = xth[4,1];
             pm[3,2] = 1 - xth[3,1];
             pm[3,3] = 1 - xth[2,1] - xth[4,1];
        endif;
    endif;
retp(pm);
endp;

/* ======================================================================= */
   @ This section echos values of parameters @

    "Order of autoregression";;pphi;
    "Order of ARCH process";;karch;
    "Number of primitive states";;ns;
    "Number of lagged states that affect y";;ps;
    "First observation used for estimation is";;nk;
    if larch == 0; "no leverage effect";  else; "with leverage effect"; endif;
    if tarch == 0;"Distribution is Normal"; elseif tarch ==1;
           "distribution is t"; endif;

proc echoo(th); @ proc to echo starting values @;
  local spar,alpha0,phi,sig0,a0,a1,g1,pm,b1,l1,cm,nu,eps,rss;

  alpha0 = th[1,1];
  spar = 2;
  if pphi > 0;
     phi = th[2:2+pphi-1,1];
     spar = spar+1;
  else;
     phi = 0;
  endif;
  "Constant term in regression";;  alpha0;
  "Autoregressive coefficients in regression";; phi';
   if barch == 0;
       sig0 = abs(th[spar,1]);
      "Initial variance";;sig0;
      spar = spar + 1;
   else;
       "Initial variance not neeeded ";
   endif;
   a0 = abs(th[spar,1]);
   a1 = abs(th[spar+1:spar+karch,1]);
  "Constant term in ARCH process";;a0;
  "Coefficients on lagged epsilon squared in ARCH process";;a1';
  spar = spar+1+karch;
  pm = matpm(th[spar:spar+(ns*(ns-1))-1-irs,1]);
  "(Transposed) matrix of transition probabilities";;pm;"";
  "The state with no adjustment to ARCH process is state 1, with transition";
  "probability ";;pm[1,1];
  spar = spar+(ns*(ns-1))-irs;

  g1 = abs(th[spar:spar+ns-2,1]);
  "Vector of variance factors for states 2 through";;ns;;g1';
  spar = spar+ns-1;

  if larch == 1;
     l1 = abs(th[spar,1]);
     "Coefficient on negative lagged change for asymmetric effect";;l1;
     spar = spar+1;
  endif;
  if tarch == 1;
       nu = 2 + abs(th[spar,1]);
       "degree of freedom for t distribution is";; nu;
       spar = spar + 1;
  endif;

retp(a0);
endp;

/* ================================================================ */
  @ This section calls main programs @

hp = pattern1;
#include procs;
x = startval;

 @ The following lines are for convenience of analysis and should be removed
             for final calculations
   izz = 2;
   kc = 2;
   ks = 2; @

call echoo(x);
"";"Initial values:";; x';
"Initial value for negative log likelihood:";; ofn(x);

"";"Do you wish to continue (y or n)?";;
  zzs = cons;
  if zzs .$== "n";
      end;
  endif;

@  goto jump; @

/* ==================================================================== */
@ Set parameters to use Gauss numerical optimizer @

     library optmum;
     #include optmum.ext;
      __btol = 1.e-08; @ This controls convergence criterion for coefficients@
      __gtol = 1.e-08; @ This controls convergence criterion for gradient @
      __algr = 1;     @ This chooses BFGS optimization  @
      __miter = 80;  @ This controls the maximum number of iterations @
      __output = 1;   @ This causes extra output to be displayed @
      __covp = 0;     @ This speeds up return from OPTMUM @

@ Next call  the GAUSS numerical optimizer @
        {x,f,g,h} =optmum(&ofn,startval);
       
"";"";"======================================================";
"          FINAL ESTIMATES";
"";"Value of log likelihood:";;-f;
"";"Coefficients:";x';"";
call echoo(x);
"";"Gradient vector:";g';
    vg=hessp(&ofn,x);
    {va,ve} = eigrs2(vg);
    va = sortc(va~ve',1);
      if va[1,1] > 0;
        "Standard errors:";
         h=invpd(vg);
         hh=sqrt(diag(h));
         hh';
       else;
          "Hessian not positive definite; eigenvalues are";
           va[.,1]';
           "eigenvector associated with smallest eigenvalue is";
           va[1,.];
       endif;

/* ======================================================================= */
@ Print out complete analysis @
jump:

kc = 2;
ks = 2;
call ofn(x);
nxx = captst; @ Use nxx = captst for full output @
"Probabilities for primitive states";

"filtered probabilities";format /rd 1,0;
"Obs ";;"return  ";;
t = 0;
do until t > ps;
  i = 1;
    do until i == ns;
       "P(st-";;t;;"=";;i;;") ";;
     i = i+1;
     endo;
  t = t+1;
endo;"";
format /rd 6,4;
 skif = (skif*hp')*(eye(ps+1).*.id[.,1:ns-1]);
 skif =  seqa(nk,1,captst)~y[nk:capt,1]~skif;
 skif[1:nxx,.];

"";"smoothed probabilities";
format /rd 1,0;
"Obs ";;"return ";;
i = 1;
   do until i > ns;
      "P(st = ";;i;;") ";;
   i = i+1;
   endo;
format /rd 6,4;
 skis = skis*hp';
 skis = seqa(nk,1,captst)~y[nk:capt,1]~skis[.,1:ns];
 skis[1:nxx,.];
"";"   Obs   Residual  Variance  Prob of observing larger value";
outprob = 1 - outprob;
iitlow = 0;
iithigh = 0;
it = 1;
do until it >  nxx;
    nk+it-1;;"   ";;resid[nk+it-1,1];;"   ";;varfor[nk+it-1,1];;"   ";;
    if outprob[it+nk-1,1] > 0.975;
            outprob[it+nk-1,1];;"  *";
            iitlow = iitlow+1;
    elseif outprob[it+nk-1,1] < 0.025;
            outprob[it+nk-1,1];;"  *";
            iithigh = iithigh+1;
    else;
            outprob[it+nk-1,1];
    endif;
it = it+1;
endo;
"";"";"Number of observations below .025 level:";iitlow;
      "Number of observations above .975 level:";iithigh;
/* hh = sortc(skis[1:captst,1]~outprob[nk:capt,1],2);
"";"Observations with lowest p values";
"";"   Obs   P-value";
format /m1 /ros 16,8;
hh[1:25,.];
"";"Observations with highest p values";  */

format /m1 /ros 16,8;

hh = skis[1:captst,1]~outprob[nk:capt,1];
hh = sortc(hh,2);
hh;