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

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

"MEXICO - Stock Returns (USD) - SWARCH(2,1)  -%";
"";
            dif = 1;        @ 1 if dep. variable to be differenced @
           intdif = 0;      @ 1 if dep. var to be a differenc of 2 vars @
           intp = 1;      @ 0 if dep. variable to be first differenced @
           mexog = 0;       @ 1 if exogeneous variable in mean @
           vexog = 0;       @ 1 if exogeneous variable in variance @
           captz = 530;     @ capt is the sample size of exog var @
           captx = 530;     @ capt is the sample size of the series@
           capt = 530;      @ capt is the sample size @
           nfil = 5;        @ number of vars. in input file @
           nfilz = 2;       @ number of vars. in input file @
            nfilx=1;        @ number of vars. in input file of exog var @


            kr = 5;           @ select the dep. var in input file @
            kz = 2;           @ select the indep. var in input file @
            kdep1 = 3;           @ select the dep. var in input file @

@     load yz[capt,nfil] = ..\data\emer\mex.2; @
@     load yz[capt,nfil] = C:\data\seb\argpes1w.prn;@
@     load yz[captz,nfilz] = ..\data\seb\arglong.prn;@
        load yz[capt,nfil] = C:\data\emer\emerd.prn;
    load yz1[captz,nfilz] = c:\data\seb\pacir.prn;

     xx = zeros(captx,1);
     zxy = xx[1:captx,1];

if intdif==1;
    y = yz[capt-captx+1:capt,kr] - yz1[.,kz];
else;
     y = yz[capt-captx+1:capt,kr];
endif;

@ yz[capt-captx+1:capt,.]~y~yz1[.,kz]; @
capt = rows(y);

/*
           x1[407,1]=1;
           x2[18,1]=1;
*/

    if dif == 1;
        y1=y;
        y=zeros(capt,1);
         if intp==1;
           yy = 100*(ln(y1[2:capt,1])-ln(y1[1:capt-1,1]));
         else;
           yy = (y1[2:capt,1] - y1[1:capt-1,1]);
         endif;
        y=yy;
        capt=capt-1;
    endif;

    if vexog>0;
           wlag=1;
           vx=abs(y[1:capt,1]);
    endif;
capt=rows(y);

@--------- Adjust any of the following to control specification desired --- @
           nk = 3;           @ nk is the first observation to be used in
                               estimation @
           ns = 2;           @ 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 = 2;          @ 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
                                  ipm = 4 user input code @
           rest = 0;         @ rest = restrictions in the prob matrix @
           karch = 0;       @ karch is the order of the ARCH process @
           karch1=1;
           larch = 0;        @ 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 = 0;        @ tarch = 0 means normal
                               tarch = 1 means t distribution
                               tarch = 2 means GED distribution @
           earch = 0;        @  especial case of EARCH -alla Pantula @
           mswarch = 0;      @  if SWARCH-M effect @
           mk = 0;          @ number of step-ahead variance forecasts @
           kvarout = 1;
           kprobout = 0;
           klb = 1;

          if ipm==2;
            rest = 0;
          endif;

          if karch>0;
          karch1=karch;
          endif;

@------ 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(2,1),...,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[6,1] =
     0.48      0.095        11.12         6.39       0.34        6.9 
;
      retp(th);
      endp;

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

ps = karch1;       @ 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 @
kau = 12;

   "Descriptive statistics of dependent variable";
capt=rows(y);
yxx=y[nk:capt,1];
captst=rows(yxx);
"Sample size:";;captst;
"mean = ";;meanc(yxx);
vyx = ((yxx-meanc(yxx))'*(yxx-meanc(yxx)));
vyxx = ((yxx-meanc(yxx))'*(yxx-meanc(yxx)))/(captst-1);
"st. deviation = ";;sqrt(vyxx);
sk = (ones(captst,1)'*(yxx-meanc(yxx))^3)/captst;
sk=sk/((vyx/captst)^1.5);
"SK= ";;sk;
ek = (ones(captst,1)'*(yxx-meanc(yxx))^4)/captst;
ek=(ek/((vyx/captst)^2))-3;
"EK= ";;ek;
nt=(captst/6)*(sk^2+.25*ek^2);
"Normality Test -Jarque-Bera (1980):";;nt;
"";
if klb==1;
"Sample autocorrelations";
"kau:";;kau;
rho=zeros(kau,1);

kx=0;
vrho=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx))/captst;

kx=1;
do until kx > kau;
rho[kx,1]=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx));
rho[kx,1]=rho[kx,1]/(captst*vrho);
kx=kx+1;
endo;
"rho:";
rho';
"";
"Standard Error:";;sqrt(1/captst);
q=0;
kq=1;
do until kq > (kau/2);
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-(kau/2):";;q*captst*(captst+2);
kq=1;
q=0;
do until kq > kau;
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-kau:";;q*captst*(captst+2);
"";
endif;

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[3,1] = 0; @
@           pm[3,3] = 0; @
           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 ns=3 @
        if izz == 1;
             pm[ns,.] = ones(1,ns);

             pm[1,1] = xth[1,1];
             pm[1,2] = xth[2,1];
             pm[1,3] = xth[3,1];
             pm[2,1] = xth[4,1];
             pm[2,2] = 1;
             pm[2,3] = xth[5,1];
             pm[3,2] = 0;

             pm = pm^2;
             pm = pm./(sumc(pm)');
        else;
        pm[ns,.] = (1 -sumc(pm))';
        endif;
     elseif ipm == 4;  @ This section can be rewritten by user to impose zeros
                          and ones where desired ns=4 @
         if izz ==1;
             pm[ns,.] = ones(1,ns);
             pm[1,1] = xth[1,1];
             pm[1,2] = xth[2,1];
             pm[1,3] = xth[3,1];
             pm[1,4] = xth[4,1];
             pm[2,1] = xth[5,1];
             pm[2,2] = xth[6,1];
             pm[2,3] = xth[7,1];
             pm[2,4] = xth[8,1];
             pm[3,1] = xth[9,1];
             pm[3,2] = 0;
             pm[3,3] = xth[10,1];
             pm[3,4] = 0;

             pm = pm^2;
             pm = pm./(sumc(pm)');

/*
             pm[1,1] = xth[1,1]^2/(1+xth[1,1]^2+xth[3,1]^2);
             pm[1,2] = xth[2,1]^2/(1+xth[2,1]^2+xth[4,1]^2);
             pm[1,3] = 0;
             pm[2,1] = xth[3,1]^2/(1+xth[1,1]^2+xth[3,1]^2);
             pm[2,2] = xth[4,1]^2/(1+xth[2,1]^2+xth[4,1]^2);
             pm[2,3] = xth[5,1]^2/(1+xth[5,1]^2);
             pm[3,1] = 1/(1+xth[1,1]^2+xth[3,1]^2);
             pm[3,2] = 1/(1+xth[2,1]^2+xth[4,1]^2);
             pm[3,3] = .5/(.5^2+xth[5,1]^2);
*/
        else;
             pm[1,1] = xth[1,1];
             pm[1,4] = 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] = xth[5,1];
             pm[3,3] = 1 - xth[4,1];
             pm[3,4] = 1 - xth[2,1];
             pm[4,2] = 1 - xth[3,1] - xth[5,1];
        endif;
/*             pm[1,1] = xth[1,1];
             pm[1,4] = xth[2,1];
             pm[2,1] = 1 - pm[1,1];
             pm[2,2] = xth[3,1];
             pm[2,3] = xth[4,1];
             pm[3,2] = xth[5,1];
             pm[3,3] = 1 - pm[2,3];
             pm[3,4] = 1 - pm[1,4];
             pm[4,2]=  1- pm[2,2] - pm[3,2];
 */
        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 karch == 0; "no ARCH effect"; else; "with ARCH effect"; endif;
    if earch == 0; "no E-SWARCH effect"; else; "with E-SWARCH effect"; endif;
    if larch == 0; "no leverage effect"; else; "with leverage effect"; endif;
    if vexog == 0; "no exog var in var"; else; "with exog var in var"; endif;
    if mswarch == 0;"No SWARCH-M"; else; "with SWARCH-M";endif;
    if tarch == 0;"Distribution is Normal"; elseif tarch ==1;
           "distribution is t"; elseif tarch ==2;
           "distribution is GED"; endif;

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

  alpha0 = th[1,1];
  "Constant term in regression ";;  alpha0;
  spar = 2;
  if mexog > 0;
     mex = th[2:2+mexog-1,1];
     spar = spar+mexog;
  "Coeffcient for Exogeneous Variable in mean ";; mex';
  else;
     mex = 0;
  endif;
  if pphi > 0;
     phi = th[spar:spar+pphi-1,1];
     spar = spar+pphi;
  else;
     phi = 0;
  endif;
  "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]);
  "Constant term in ARCH process";;a0;
if karch >0;
   a1 = abs(th[spar+1:spar+karch,1]);
  "Coefficients on lagged epsilon squared in ARCH process";;a1';
endif;
  spar = spar+1+karch;

  pm = matpm(th[spar:spar+(ns*(ns-1))-1-rest,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))-rest;
  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;
    if earch == 1;
     l1 = th[spar,1];
    else;
     l1 = abs(th[spar,1]);
    endif;
     "Coefficient on negative lagged change for asymmetric effect";;l1;
     spar = spar+1;
  endif;

  if vexog > 0;
     vex = abs(th[spar:spar+vexog-1,1]);
  "Coeffcient for Exogeneous Variable in variance ";; vex';
  else;
     vex = 0;
  endif;
     spar = spar + vexog;

  if mswarch > 0;
     msw = abs(th[spar:spar+mswarch-1,1]);
  "Coeffcient for SWARCH in mean ";; msw';
  else;
     msw = 0;
  endif;
     spar = spar + mswarch;

  if tarch == 1;
       nu = 2 + abs(th[spar,1]);
       "degree of freedom for t distribution is";; nu;
       spar = spar + 1;
  endif;
  if tarch == 2;
       nu =abs(th[spar,1]);
       "GED distribution parameter is";; nu;
       spar = spar + 1;
  endif;
retp(a0);
endp;

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

hp = pattern1;
#include c:\prog\eprocs1.txt;
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;

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

 
     library optmum;
     #include optmum.ext;
      _opbtol = 1.e-06; @ This controls convergence criterion for coefficients@
      _opgtol = 1.e-06; @ This controls convergence criterion for gradient @
      _opalgr = 1;     @ This chooses BFGS optimization  @
      _opmiter = 50;  @ This controls the maximum number of iterations @
      _opstmth = "steep stepbt";
      _opmdmth = "bfgs golden";
   /*   __covp = 0;  */   @ This speeds up return from OPTMUM @

@ Next call  the GAUSS numerical optimizer @
        output file=c:\emer\junk.out off;
        output file=junk reset;
        {x,f,g,h} =optmum(&ofn,startval);
       output file=junk off;
        output file=c:\emer\junk.out on;

"";"";"======================================================";
"          FINAL ESTIMATES";
"";"Value of log likelihood:";;-f;
"";
call echoo(x);
"";"Coefficients:";x';"";
"Gradient vector:";g';
@ goto jump; @
    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';
        "T-statistics:";
         tst=x./hh;
         tst';
       else;
          "Hessian not positive definite; eigenvalues are";
           va[.,1]';
           "eigenvector associated with smallest eigenvalue is";
           va[1,.];
        endif;

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


kc = 2;
ks = 2;
call ofn(x);

nxx = captst; @ Use nxx = captst for full output @

if kvarout == 1;
output file=c:\emer\junk.out off;
@ output file=junk.fil reset; @
/*
"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,.];
output file=junk.fil off;
*/
output file=c:\emer\junk.smo reset;
format /rd 6,4;
 skis = skis*hp';
 skis = seqa(nk,1,captst)~y[nk:capt,1]~skis[.,1:ns];
 skis[1:nxx,.];
output file=c:\emer\junk.smo off;

output file=junk.var reset;
"";"   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;
output file=junk.var off;

output file=c:\emer\junk.out on;
"";"";"Number of observations below .025 level:";iitlow;
      "Number of observations above .975 level:";iithigh;
endif;

/*
"";
#include nnproc;
yxx=resid[nk:capt,1]./sqrt(varfor[nk:capt,1]);
xxx=zeros(capt-pphi,1);

    iz = 1;
    do until iz > pphi;
       xxx=xxx~y[pphi+1-iz:capt-iz,1];
    iz = iz+1;
    endo;

neural(yxx,xxx[nk-pphi:capt-pphi,2:pphi+1],pphi,0);
*/
"";
msx=resid[nk:nxx,1]^2-varfor[nk:nxx,1];
"MSE:";;sumc(msx^2)/nxx;
"MAE:";;sumc(abs(msx))/nxx;
lsx=ln(resid[nk:nxx,1]^2)-ln(varfor[nk:nxx,1]);
"ALES:";;sumc(lsx^2)/nxx;
"ALE:";;sumc(abs(lsx))/nxx;

/* 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";  */

output file=junk.out off;

if kprobout ==1;
format /m1 /ros 16,8;
output file=junk.prob reset;
hh = skis[1:captst,1]~outprob[nk:capt,1];
hh = sortc(hh,2);
hh;
output file=junk.prob off;
endif;

output off;