/* MARCH  */
/*         GENERAL ARCH PROGRAM
           This is a Gauss program.  It  chooses parameters so as to maximize
           the likelihood function for a STAR in mean with ARCH.
           (the negative of the log likelihood is calculated by a separate
           Gauss procedure, contained in the file GLARCH).
*/

                   output file=tango reset;

     library optmum;
     #include optmum.ext;

      _opbtol = 1.e-15; @ This controls convergence criterion for coefficients@
      _opgtol = 1.e-15; @ This controls convergence criterion for gradient @
      _opalgr = 2;@  This chooses Davidon-Fletcher-Powell optimization  @
      _opstmth = "steep one";
      _opmdmth = "bfgs steptb";
      _opmiter = 105;  @ This controls the maximum number of iterations @

"-------------------------- OPTIMIZATION RESULTS  --------------------------";

/*
    Variable meanings:
            dif = 0 if you want to log difference the data
            n = number of observations
            nfil = number of files in load matrix
            kr = file selected
            nk = number of parameters
            y = data to be analyzed
            exog = if exogenous variables in mean included
            kexog = number of exogenous variables in mean
            barch = 0 if initial observation estimated by MLE
            karch = arch lags
            garchb= garch terms
            larch = leverage effect
            vexog = if exogenous variables in variance included
            kvexog = number of exogenous variables in variance
            ki = initial observation to start iteration
            earch = exponential garch formulation
            kc = 1 suppresses output print
            kd = distribution (1=normal, 2=ged, 3=student)
            km = 1;   1=garch(1,1)  2=arch(2)  3=ding  4=egarch
   Don't forget to change the parameter dimension in STARTVAL below by hand */

    format /rd 15,6;
 /* Specification of the model */
    nx = 996;      @ sample size @
    n1 = 996;
    n = 996;
    nfil = 3;     @ number of files in  data file @
    nfilx = 1;     @ number of files in exog. data file @
    nfilz = 3;     @ number of files in exog.  data file @
    kr = 3;       @ column for dependent variable @
    kr2 = 1;
    kz = 2;
    dif = 0;      @ 1 if differenced data @
    earch = 0;    @ 1 if EGARCH model @
    larch = 0;    @ 1 if leverage effect @
    pphi = 0;
    exog = 0;     @ 1 if exogenous variable in mean @
    vexog = 0;    @ 1 if exogenous variable in variance @
    kexog = 0;
    kvexog = 0;
    mnl = 0;      @ 1 if STAR effect in mean @
    barch = 2;
    karch = 1;
    garchb = 1;
    nlarch = 0;   @ if STAR in variance @
    nkarch = 0;
    kmean = 0;
    kvarp = 1;  @=1 if printing var in tango.var file @
    ksterr = 1;  @=1 if printing standard errors in tango output file @
    mk = 0;     @ nth step ahead out-of-sample forecast @
    ki = 4;     @ initial observation to be used @
    kd = 1;      /* 1=normal   2=ged   3=student */
    km = 1;
    klb = 0;    @ print lbj for raw series @
    kau =10;
    if earch ==1;
    km=3;
    endif;
    captst = n1-ki+1;

"GBP/USD-  E-GARCH(1,1) * GED"
"";
   "Descriptive statistics of dependent variable";
     load yz[nx,nfil] = data1.mc;
     load zx[nx,nfilx] = ..\data\ret\world.ret;
     load yyi[n1,nfilz] = ..\data\ret\eur.w;
     yi = 100*ln(1+yyi[.,kz]*.01)/52;
      y = yz[.,kr];
     xx=zx[nx-n1+1:nx,kr2]-yi;

     proc startval;
        local th;

  let th[4,1] =
  0.01255   0.00785   0.07644   0.89746
;
         retp(th);
           endp;

if dif == 1;
     yy = (ln(y[2:n1,1])-ln(y[1:n1-1,1]))*100;
     n1=n1-1;
     captst = n1-ki+1;
     y=zeros(n1,1);
     y=yy;
endif;
"captst:";;captst;

if vexog == 1;
    xv=abs(y);
endif;

yxx=y[ki:n,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;
"*Level series";
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);
"";
"*Squared series";
rho=zeros(kau,1);
yxx=yxx^2;
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';
"";
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;
/* =================================================================== */

kc = 1;

/* This section echos values of parameters */

    "Order of autoregression ";;pphi;
    if exog ==1; "Exogenous variable in mean ";endif;
    if vexog ==1; "Exogenous variable in variance ";endif;
    "Order of ARCH process ";;karch;
    "Order of GARCH process ";;garchb;
    "First observation used for estimation is ";;ki;
    if mnl == 0; "no STAR in mean"; else; "with mean STAR effect"; endif;
    if earch == 0; "no exponential garch";  else; "EGARCH model"; endif;
    if larch == 0; "no leverage effect";  else; "with leverage effect"; endif;
    if nlarch == 0; "no NLARCH";  else; "with NLARCH effect"; endif;
    if kd == 1;"Distribution is Normal"; elseif kd ==3;
           "distribution is t"; elseif kd ==2;
           "distribution is GED"; endif;

proc echoo(th); @ proc to echo starting values @;
  local spar,alpha0,phi,beta,sig0,a0,a1,g1,pm,b1,l1,vbeta,cm,nu,eps,
    rss,nlar0,nlar1,malpha0,mphi,mlm,mlc,nlc,nlm,pex;

  alpha0 = th[1,1];
  "Constant term in regression";;  alpha0;
  spar = 2;
  if pphi > 0;
     phi = th[2:2+pphi-1,1];
     spar = spar+pphi;
  "Autoregressive coefficients in regression";; phi';
  else;
     phi = 0;
  endif;
    if exog >0;
    beta = th[spar:spar+exog-1,1];
  "Exogenous variable coefficients in mean";; beta';
     spar = spar+exog;
    endif;
   if mnl == 1;
    malpha0 = th[spar,1];
  "Constant NL term in regression";;  malpha0;
    spar = spar+1;
   if pphi > 0;
     mphi = th[spar:spar+pphi-1,1];
     spar = spar+pphi;
  "Autoregressive NL coefficients in regression";; mphi';
  else;
     phi = 0;
  endif;
     mlc = th[spar,1];
     "Coefficient on constant in logistic-mean part";;mlc;
     spar = spar+1;
     mlm = th[spar,1];
     "Coefficient on mean in logistic-mean part";;mlm;
     spar = spar + 1;
   endif;

   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;

  if garchb >0;
    b1=th[spar,1];
  "Coefficients on lagged variance in GARCH process";;b1';
    spar=spar+garchb;
  endif;

  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 km==4;
     pex = abs(th[spar,1]);
     "Exponenet coefficient on Ding et al. model:";;pex;
      spar = spar + 1;
    endif;

    if vexog ==1;
    vbeta = th[spar:spar+kvexog-1,1];
  "Exogenous variable coefficients in variance";; vbeta';
     spar = spar+kvexog;
    endif;

    if kmean ==1;
    cm = th[spar:spar+1-1,1];
  "ARCH in mean coefficients ";; cm';
     spar = spar+1;
    endif;

 if nlarch == 1;
     nlar0 = abs(th[spar,1]);
     nlar1 = abs(th[spar+1:spar+nkarch,1]);
     "Coefficient on constant in NLARCH part";;nlar0;
     "Coefficient on lagged ARCH effect in NLARCH part";;nlar1';
     spar = spar+1+nkarch;
     nlc = th[spar,1];
     "Coefficient on constant in logistic-NLARCH part";;nlc;
     spar = spar+1;
     nlm = th[spar,1];
     "Coefficient on mean in logistic-NLARCH part";;nlm;
     spar = spar + 1;
 endif;

 if kd == 3;
       nu = 2 + abs(th[spar,1]);
       "degree of freedom for t distribution is";; nu;
       spar = spar + 1;
  endif;

  if kd == 2;
       nu =abs(th[spar,1]);
       "GED distribution parameter is";; nu;
       spar = spar + 1;
  endif;
retp(a0);
endp;

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

        "";"";


#include larch;
x = startval;

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;

        output off;
        {x,f,g,retcode} =optmum(&ofn,startval);
        output file=tango on;


"";"";"          FINAL ESTIMATES";
"";"Value of negative log likelihood: ";;f;
call echoo(x);
"";"Final parameters: ";x';
"";"Gradient vector: ";g';
if ksterr==1;
    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';
       endif;
endif;
"";
        kc = 2;
        call ofn(x);
        kc =1;
        output off;