@ This GAUSS file reads in data, sets options, and calls numerical
    optimization routine for numerical estimation of a Bivariate SWARCH model  with state 
    dependent correlations. Reference Ramchand and Susmel (1998) J. of Empirical Finance. @

        output file=junkw.var reset;
        output file=junks.out reset;
        format /m1 /ros 11,5;
/* ======================================================================= */
@  Alter this section so as to read in your data and set parameters @
"HONG - ARG Stock Market Weekly Multi-SWARCH(2,1,2) ";
"";
           @ Crisp data @
            dif1 = 0;
            dif2 = 1;
            ki = 504;
            ki2 = 0;
            nm=2;             @ number of series in multivariate model @
           capt = 1033;        @ capt is the sample size @
           capt2 = 530;        @ capt is the sample size @
          load yy1[capt,1] = c:\data\ret\hong.ret;
@           load yy1[capt2,5] = c:\data\emer\emerd.prn;@
           load yy2[capt2,5] = c:\data\emer\emerd.prn;
           x1=yy1[ki+1:capt,1];
           x2=yy2[ki2+1:capt2,2];

y1=x1;
y2=x2;
capt = rows(y1);
capt2 = rows(y2);

    if dif1 == 1;
        y=y1;
        y1=zeros(capt-1,1);
         yy1=100*(ln(y[2:capt,1])-ln(y[1:capt-1,1]));
        y1=yy1;
        capt=capt-1;
    endif;
capt=rows(y1);

    if dif2 == 1;
        y=y2;
        y2=zeros(capt2-1,1);
         yy2=100*(ln(y[2:capt2,1])-ln(y[1:capt2-1,1]));
        y2=yy2;
        capt2=capt2-1;
    endif;

           x1 = zeros(capt,1);
           x2 = zeros(capt,1);
           x1[90,1]=1;
           x2[107,1]=1;
           xx=x1~x2;

@--------- Adjust any of the following to control specification desired --- @
           nk = 4;           @ nk is the first observation to be used in
                               estimation @
           mexog1 = 0;        @ 1 if exogeneous variable in mean 1 @
           mexog2 = 0;
           ns1 = 2;
           ns2 = 2;
           ns = 4;           @ ns is the number of primitive states @
           pphi1 = 1;        @ pphi is the number of lags in autoregression
                                for y; @
           pphi2 = 1;
           izz = 1;          @ izz = 1 means pij parameterized as
                                       th(ij)^2/sum(th(ij)^2)
                               izz = 2 means pij parameterized as pij @
           ipm = 5;          @ 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 = 6;         @ rest = restrictions in the prob matrix @
           karch1 = 1;       @ karch1 is the order of the ARCH process 1 @
           kzarch1 = 1;      @ kzarch1 is the order of the ARCH process 1 @
           karch2 = 1;       @ karch2 is the order of the ARCH process 2 @
           kzarch2 = 1;      @ kzarch2 is the order of the ARCH process 2 @
           larch1 = 0;       @ larch = 0 means no leverage effect,
                               larch = 1 means leverage effect @
           larch2 = 0;
           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;        @  EARCH @
           cov = 1;          @ if covariance calculated @
           lam =3;
           kvar = 0;

        if karch1>0;
        kzarch1=karch1;
        endif;

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

        if ipm==4;
        rest=8;
        endif;

        if ipm==5;
        rest=10;
        endif;

        if ipm==7;
        rest=6;
        endif;


        if karch2>0;
        kzarch2=karch2;
        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[25,1]=
 0.47915   0.40338 -0.038261  0.024767    6.8539 2.7114E-06   -10.164 
  0.24097    451.43   0.91888 0.00049929   0.30009   -41.206    7.0617 
0.0051840   0.27320   -56.336    .0001 2.7914  -0.42194    5.3186    10.684 
 2.134
;
      retp(th);
      endp;

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

ps = kzarch1;  @ ps is the number of lagged states that matter @
n = ns^(ps+1);     @ n is the dimension of the state vector @
kk=(n/2);
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,nm);  @ 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 @
"";

   "Descriptive statistics of dependent variable 1";
yxx=y1[nk:capt,1];
"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;
"";
   "Descriptive statistics of dependent variable 2";
yxx=y2[nk:capt,1];
"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;
"";

/* Useful matrix */
    hw1=zeros(4,8);
    hw2=zeros(4,8);
    hh=zeros(8,8);
    hw1[.,1]=1|0|0|0;
    hw1[.,2]=0|0|0|0;
    hw1[.,3]=1|0|0|0;
    hw1[.,4]=0|0|0|0;
    hw1[.,5]=0|1|0|0;
    hw1[.,6]=0|0|0|0;
    hw1[.,7]=0|1|0|0;
    hw1[.,8]=0|0|0|0;

    hw2[.,1]=0|0|0|0;
    hw2[.,2]=1|0|0|0;
    hw2[.,3]=0|0|0|0;
    hw2[.,4]=0|1|0|0;
    hw2[.,5]=0|0|0|0;
    hw2[.,6]=1|0|0|0;
    hw2[.,7]=0|0|0|0;
    hw2[.,8]=0|1|0|0;

hww1=zeros(2,8)|hw1;
hww2=zeros(2,8)|hw2;

    hh[.,1]=1|0|1|0|1|0|1|0;
    hh[.,2]=0|1|0|1|0|1|0|1;
    hh[.,3]=1|0|1|0|1|0|1|0;
    hh[.,4]=0|1|0|1|0|1|0|1;
    hh[.,5]=1|0|1|0|1|0|1|0;
    hh[.,6]=0|1|0|1|0|1|0|1;
    hh[.,7]=1|0|1|0|1|0|1|0;
    hh[.,8]=0|1|0|1|0|1|0|1;

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,px11,px12,px21,px22,ixth,py11,py12,py21,py22;
   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[1,4] = 0;
             pm[2,1] = xth[4,1];
             pm[2,2] = 1;
             pm[2,3] = 0;
             pm[2,4] = xth[5,1];
             pm[3,1] = 1;
             pm[3,2] = 0;
             pm[3,3] = xth[6,1];
             pm[3,4] = 1;
             pm[4,1] = 0;
             pm[4,2] = 0;
             pm[4,4] = 0;
             @ pm[2,3] = xth[5,1]; @
             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;
             px11= xth[1,1]^2/(1+xth[1,1]^2);
             px12 = (1-px11);
             px21= xth[2,1]^2/(1+xth[2,1]^2);
             px22= (1-px21);
             py11= xth[3,1]^2/(1+xth[3,1]^2);
             py12 = (1-py11);
             py21= xth[4,1]^2/(1+xth[4,1]^2);
             py22= (1-py21);

             pm[1,1] =px11*py11;
             pm[1,2] =px11*py21;
             pm[1,3] =px21*py11;
             pm[1,4] =px21*py21;
             pm[2,1] =px11*py12;
             pm[2,2] =px11*py22;
             pm[2,3] =px21*py12;
             pm[2,4] =px21*py22;
             pm[3,1] =px12*py11;
             pm[3,2] =px12*py21;
             pm[3,3] =px22*py11;
             pm[3,4] =px22*py21;
             pm[4,1] =1-pm[1,1]-pm[2,1]-pm[3,1];
             pm[4,2] =1-pm[1,2]-pm[2,2]-pm[3,2];
             pm[4,3] =1-pm[1,3]-pm[2,3]-pm[3,3];
             pm[4,4] =1-pm[1,4]-pm[2,4]-pm[3,4];
        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;
     elseif ipm == 5;
         if izz ==1;
             px11= xth[1,1]^2/(1+xth[1,1]^2);
             px12 = (1-px11);
             px21= xth[2,1]^2/(1+xth[2,1]^2);
             px22= (1-px21);

             pm[1,1] =px11;
             pm[1,2] =0;
             pm[1,3] =0;
             pm[1,4] =px21;
             pm[2,1] =0;
             pm[2,2] =0;
             pm[2,3] =0;
             pm[2,4] =0;
             pm[3,1] =0;
             pm[3,2] =0;
             pm[3,3] =0;
             pm[3,4] =0;
             pm[4,1] =px12;
             pm[4,2] =0;
             pm[4,3] =0;
             pm[4,4] =px22;
        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;
     elseif ipm == 6;
         if izz ==1;
             px11= xth[1,1]^2/(1+xth[1,1]^2);
             px12 = (1-px11);
             px21= xth[2,1]^2/(1+xth[2,1]^2);
             px22= (1-px21);
             py11= xth[3,1]^2/(1+xth[3,1]^2);
             py12 = (1-py11);
             py21= xth[4,1]^2/(1+xth[4,1]^2);
             py22= (1-py21);

             pm[1,1] =px11*py11;
             pm[1,2] =px11*py21;
             pm[1,3] =px21*py11;
             pm[1,4] =px21*py21;
             pm[2,1] =px11*py12;
             pm[2,2] =px11*py22;
             pm[2,3] =px21*py12;
             pm[2,4] =px21*py22;
             pm[3,1] =px12*py12;
             pm[3,2] =px12*py22;
             pm[3,3] =px22*py22;
             pm[3,4] =px22*py22;
             pm[4,1] =1-pm[1,1]-pm[2,1]-pm[3,1];
             pm[4,2] =1-pm[1,2]-pm[2,2]-pm[3,2];
             pm[4,3] =1-pm[1,3]-pm[2,3]-pm[3,3];
             pm[4,4] =1-pm[1,4]-pm[2,4]-pm[3,4];
        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;
     elseif ipm == 7;
         if izz ==1;
             pm[ns,.] = ones(1,ns);
             pm[1,1] = xth[1,1];
             pm[1,2] = xth[2,1];
             pm[1,3] = 0;
             pm[1,4] = xth[3,1];
             pm[2,1] = xth[4,1];
             pm[2,2] = xth[5,1];
             pm[2,3] = 0;
             pm[2,4] = xth[6,1];
             pm[3,1] = 0;
             pm[3,2] = 0;
             pm[3,3] = .9;
             pm[3,4] = 0;
             pm[4,3] = .1;
             pm = pm^2;
             pm = pm./(sumc(pm)');

        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;

        endif;
retp(pm);
endp;

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

    "Order of autoregressions";;pphi1~pphi2;
    "Order of ARCH processes";;karch1~karch2;
    "Number of primitive states";;ns;
    "Number of lagged states that affect y";;ps;
    "First observation used for estimation is";;nk;
   if earch == 0; "no E-SWARCH effect"; else; "with E-SWARCH effect"; endif;
   if larch1 == 0; "no leverage effect"; else; "with leverage effect"; endif;
   if tarch == 0;"Distribution is Normal"; elseif tarch ==1;
           "distribution is t"; elseif tarch ==2;
           "distribution is GED"; endif;
   if ipm == 4; "independent states model"; else; "no restrictions on P"; endif;
"";
proc echoo(th); @ proc to echo starting values @;
  local spar,alpha0,alpha02,mex1,mex2,phi,sig0,a0,a1,g1,pm,b1,l1,cm,nu,eps,
  rss,a01,a02,a11,a12,phi1,phi2,corr,l2,g2;

  alpha0 = th[1,1];
  alpha02 = th[2,1];
  "Constant term in regression ";; alpha0~alpha02;
  spar = 3;
  if mexog1 > 0;
     mex1 = th[spar:spar+mexog1-1,1];
     spar = spar+mexog1;
  "Coeffcient for Exogeneous Variable in mean 1";; mex1;
  else;
     mex2 = 0;
  endif;

  if mexog2 > 0;
     mex2 = th[spar:spar+mexog2-1,1];
     spar = spar+mexog2;
  "Coeffcient for Exogeneous Variable in mean 2";; mex2;
  else;
     mex2 = 0;
  endif;

  if pphi1 > 0;
     phi1 = th[spar:spar+pphi1-1,1];
     spar = spar+pphi1;
  else;
     phi1 = 0;
  endif;
  "Autoregressive coefficients 1 in regression";; phi1';
  if pphi2 > 0;
     phi2 = th[spar:spar+pphi2-1,1];
     spar = spar+pphi2;
  else;
     phi2 = 0;
  endif;
  "Autoregressive coefficients 2 in regression";; phi2';
   if barch == 0;
       sig0 = abs(th[spar,1]);
      "Initial variance";;sig0;
      spar = spar + 1;
   else;
       "Initial variance not neeeded ";
   endif;

   a01 = abs(th[spar,1]);
  "Constant term in ARCH process 1";;a01;
   spar = spar+1;

   if karch1 >0;
   a11 = abs(th[spar:spar+karch1-1,1]);
  "Coefficients on lagged epsilon squared in ARCH process 1:";;a11';
   spar = spar+karch1;
   endif;

   a02 = abs(th[spar,1]);
  "Constant term in ARCH process 2";;a02;
   spar = spar+1;

    if karch2>0;
   a12 = abs(th[spar:spar+karch2-1,1]);
  "Coefficients on lagged epsilon squared in ARCH process 2:";;a12';
  spar = spar+karch2;
   endif;

  pm = matpm(th[spar:spar+1*(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;
  "spar:";;spar;
  g1 = abs(th[spar:spar+ns1-2,1]);
  "Vector of variance factors for 1 states 2 through";;ns1;;g1';
  spar = spar+ns1-1;

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

  if larch1 == 1;
    if earch == 1;
     l1 = th[spar,1];
    else;
     l1 = abs(th[spar,1]);
    endif;
     "Coefficient on negative lagged change for asymmetric effect 1:";;l1;
     spar = spar+1;
  endif;

  if larch2 == 1;
    if earch == 1;
     l2 = th[spar,1];
    else;
     l2 = abs(th[spar,1]);
    endif;
     "Coefficient on negative lagged change for asymmetric effect 2:";;l2;
     spar = spar+1;
  endif;
  if cov > 0;
@    corr= th[spar:spar+cov-1,1];@
 corr = 1/(1+th[spar:spar+cov-1,1]^2);
     "Correlation coefficient:";;corr;
    spar=spar+cov;
  endif;

  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(a01);
endp;

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

hp = pattern1;
ns=ns1;
n = ns1^(ps+1);
hp1 = pattern1;
ns=ns2;
n = ns2^(ps+1);
hp2 = pattern1;
ns=4;
n = ns^(ps+1);
#include eprmul2.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-10; @ This controls convergence criterion for coefficients@
      _opgtol = 1.e-10; @ This controls convergence criterion for gradient @
      _opalgr = 1;     @ This chooses BFGS optimization  @
      _opmiter = 90;  @ 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=junks.out off;
        output file=junk reset;
        {x,f,g,h} =optmum(&ofn,startval);
        output file=junk off;
        output file=junks.out on;

"";"";"======================================================";
"          FINAL ESTIMATES";
"";"Value of log likelihood:";;-f;
"";
call echoo(x);
"";"Coefficients:";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';
        "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);

output off;