/* LARCH

   This proc evaluates the STAR likelihood function with
   an ARCH specification.
   You can select an ARCH(1), ARCH(2) or ARCH(1)-M specification.
   This proc can allow for an asymmetric effect of negative returns in
   variance -Glosten et al. (1989).
   */

/* ------------------------------------------------------- */

proc ofn(th);
 local th,alpha0,sig0,f,qq,it,phi,beta,h0,a0,a1,a2,b1,l1,vbeta,z0,z1,h1,
 spar,eta2,heta,iz,itl,ith,hsig,nheta,nhsig,nlar0,nlar1,nlc,nlm,nf,
 malpha0,mphi,mlm,mlc,mnz,msig,mnf,chqq,roundqq,k,q,kq,rho,vrho,pex,
 junk,fj,lev,cm,pm,nu,gp,lambda,lambada,x,yy,prob,v0,lx,ly,a3,sab,hf,j,
 z00,vyx,vyxx,yxx,ek,sk,nt,msx,lsx,mzz,serr,aerr,rsq,ess,aic,cged1,cnor;
   junk = zeros(16,2);
     /*  Establish values of constant parameters */

@ -------------- Convert parameter vector to convenient form --------------@
  spar = 1;
  alpha0 = th[1,1];
  spar = spar + 1;
  if pphi > 0;
     phi = th[spar:spar+pphi-1,1];
  else;
     phi = 0;
  endif;
  spar = spar + pphi;

  if exog>0;
    beta=th[spar:spar+exog-1,1];
  else;
    beta=0;
  endif;
   spar =spar+exog;

  if mnl==1;
    malpha0 = th[spar,1];
    mphi = th[spar+1:spar+pphi,1];
    spar = spar + pphi +1;
    mlc = abs(th[spar,1]);
    spar = spar + 1;
    mlm = th[spar,1];
    spar = spar + 1;
  endif;

  if barch == 0;
      sig0 = sqrt(abs(th[spar,1]));
      spar = spar + 1;
  else;
      sig0 = 10000;   @ this should not be used with this option.  The absurd
                          large value should flag if used inadvertently @
  endif;

 if earch==1;
  a0 = th[spar,1];
 else;
  a0 = abs(th[spar,1]);
 endif;
 spar=spar+1;

 if karch > 0;
 if earch==1;
  a1 = th[spar:spar+karch-1,1];
 else;
  a1 = abs(th[spar:spar+karch-1,1]);
 endif;
  spar = spar+karch;
 endif;

  if garchb==1;
   if earch==1;
        b1 = th[spar,1];
    else;
        b1 = abs(th[spar,1]);
   endif;
  else;
    b1 = 0;
    spar = spar+1;
 endif;

    l1=0;
  if larch == 1;
   if earch == 1;
     l1 = th[spar,1];
   else;
     l1 = abs(th[spar,1]);
   endif;
     spar = spar+1;
  endif;

   if km==4;
    pex = abs(th[spar,1]);
    spar = spar+1;
    else;
    pex = 2;
   endif;

  if vexog==1;
   if earch==1;
    vbeta=th[spar:spar+kvexog-1,1];
   else;
    vbeta=abs(th[spar:spar+kvexog-1,1]);
   endif;
     spar =spar+kvexog;
  endif;

  if kmean==1;
    cm=th[spar:spar,1];
     spar =spar+1;
    else;
    cm=0;
  endif;

  if nlarch==1;
    nlar0 = th[spar,1];
    nlar1 = th[spar+1:spar+nkarch,1];
    spar = spar + nkarch +1;
    nlc = abs(th[spar,1])*100;
    spar = spar + 1;
    nlm = th[spar,1];
    spar = spar + 1;
  endif;

  if kd == 2;
      gp = abs(th[spar,1]);
    lambda = sqrt(gamma(1/gp)*2^(-2/gp)/gamma(3/gp));
    lambada = gp/(lambda*gamma(1/gp)*(2^(1/gp))*2);
       spar = spar + 1;
  endif;

  if kd == 3;
      nu = 2 + abs(th[spar,1]);
    lambada = (1/sqrt(((nu - 2)*pi)))*gamma(abs((nu+1)/2))/gamma(abs(nu/2));
    spar = spar + 1;
  endif;

@ ---------- Convert data to AR resids ----------@
  z0 = y[pphi+1:n,1] - alpha0 ;

  if exog >0;
  z0 = z0 - xx[pphi+1:n,.]*beta;
  endif;

    iz = 1;
    do until iz > pphi;
       z0 = z0 - phi[iz,1]*y[pphi+1-iz:n-iz,1];
    iz = iz+1;
    endo;
    if barch == 1;
    sig0=sqrt((z0[ki:n-pphi,1]'*z0[ki:n-pphi,1])/(captst-pphi-1));
    endif;
    if pphi > 0;
       z0 = (sig0*ones(pphi,1)) | z0;
    endif;

 if mnl ==1;
    mnz = zeros(n,1);
    mnz[1+pphi:n,1] = malpha0*ones(n-pphi,1);
      mnz[1+pphi:n,1] = mnz[1+pphi:n,1] + mphi[1,1]*xx[pphi+1:n,1];

    if pphi > 0;
       mnz = (sig0*ones(pphi,1)) | mnz[1+pphi:n,1];
    endif;

    mnf=ones(n-1,1)./
    (1*ones(n-1,1) + exp(-mlc*(xx[1:n-1,1] - mlm*ones(n-1,1))));
    mnf=1/(1+ exp(-mlc*(sig0- mlm)))|mnf;
    mnz=mnz.*mnf;

    z0 = z0 + mnz;
endif;

/*
 if mnl ==1;
    mnz = zeros(n,1);
    mnz[1+pphi:n,1] = malpha0*ones(n-pphi,1);
      iz = 1;
      do until iz > pphi;
      mnz[1+pphi:n,1] = mnz[1+pphi:n,1] + mphi[iz,1]*z0[pphi+1-iz:n-iz,1];
      iz = iz + 1;
      endo;

    if pphi > 0;
       mnz = (sig0*ones(pphi,1)) | mnz[1+pphi:n,1];
    endif;

    mnf=ones(n-1,1)./
    (1*ones(n-1,1) + exp(-mlc*(z0[1:n-1,1] - mlm*ones(n-1,1))));
    mnf=1/(1+ exp(-mlc*(sig0- mlm)))|mnf;
    mnz=mnz.*mnf;

    z0 = z0 + mnz;
endif;
*/

@------------- Collect karch+1 most recent squared residuals in
               the captst by karch + 1 matrix heta ------@
  if km <4;

    eta2 = z0^2;
    heta = a0*ones(captst,1);
    h0=zeros(n,1);
  if barch == 0;
    h0[1:ki,1]=sig0^2*ones(ki,1);
  else;
    h0[1:ki,1]=eta2[1:ki,1];
  endif;

      iz = 1;
      do until iz > karch;
        if earch == 1;
      heta = heta ~ abs(z0[ki-iz:n-iz,1]);
        else;
      heta = heta ~ eta2[ki-iz:n-iz,1];
        endif;
      iz = iz + 1;
      endo;

     if larch == 1;  @ Correct heta for leverage effects if desired @
     if earch ==1;
     heta[.,2] = heta[.,2]
                + (l1/a1[1,1])* z0[ki-1:n-1,1];
    else;
@     heta[.,2] = heta[.,2] + (l1/a1[1,1])*((z0[ki-1:n-1,1] .< 0).* @
     heta[.,2] = heta[.,2] + (l1/a1[1,1])*((z0[ki-1:n-1,1] .> 0).*
                   eta2[ki-1:n-1,1]);
     endif;
     endif;
 endif;

 if km==4;

    eta2 = (abs(z0) - l1*z0)^pex;
    heta = a0*ones(captst,1);
    h0=zeros(n,1);
  if barch == 0;
    h0[1:ki,1]=sig0^pex*ones(ki,1);
  else;
    h0[1:ki,1]=eta2[1:ki,1];
  endif;

      iz = 1;
      do until iz > karch;
      heta = heta ~ eta2[ki-iz:n-iz,1];
      iz = iz + 1;
      endo;
endif;

if karch > 0;
    hsig=1|a1;
    if earch ==1;
    heta=heta*hsig-a0;
    else;
    heta=heta*hsig;
    endif;

 if nlarch ==1;
    nheta = nlar0*ones(captst,1);

      iz = 1;
      do until iz > karch;
         nheta = nheta ~ eta2[ki-iz:n-iz,1];
      iz = iz + 1;
      endo;

    nhsig=1|nlar1;
    nheta=nheta*nhsig;
       nf=ones(captst,1)./
    (1*ones(captst,1) + exp(-nlc*(abs(z0[ki-1:n-1,1]^2)-nlm*ones(captst,1))));
    nheta=nheta.*nf;
     heta = heta + nheta;
 endif;

endif;

    if vexog==1;
     heta = heta + xv[ki-1:n-1,kvexog]*vbeta;
    endif;

   /* Basic iteration */

/* GARCH(karch,1)-M with leverage GJR effect */

f = 0;
prob=zeros(n,1);
v0=zeros(n,1);
it = ki;

if kd==2;
cged1= a0 - a1*(lambda*2.0^(1/gp)*gamma(2/gp)/gamma(1/gp));
else;
cged1=a0 - a1*sqrt(2/pi);
endif;

if garchb==1;
do until it > n;

if km==1;
 h0[it,1] = heta[it-ki+1,1] + b1*h0[it-1,1];
endif;

if km==3;
 h0[it,1] = exp( heta[it-ki+1,1]/sqrt(h0[it-1,1])
 + cged1 + b1*ln(h0[it-1,1]) );
endif;


if km==4;
 h0[it,1] = (heta[it-ki+1,1] + b1*(h0[it-1,1])^pex)^1/pex;
endif;

it = it+1;
endo;
else;
h0[ki:n,1]=heta[1:captst,1];
endif;

if karch==0;
h0[ki:n]=heta[it-ki+1:captst];
endif;

if kmean==1;
z0 = z0 - cm*h0;
eta2=z0^2;
endif;


 /* Likelihood */
if kd==1;
        qq= abs(exp(-(eta2[ki:n,1])./(2*h0[ki:n,1]))./sqrt(h0[ki:n,1]));

   roundqq=.01E-120*ones(captst,1);
    chqq = maxc(qq'|roundqq');
        f = sumc(ln(chqq)) - captst*.9189385;
@        f = sumc(ln(qq)) - captst*.9189385; @
endif;

if kd==2;
    qq = exp(-.5*(abs(z0[ki:n,1]./(sqrt(h0[ki:n,1])*lambda)))^gp)
    ./sqrt(h0[ki:n,1]);
    roundqq=.01E-70*ones(captst,1);
    chqq = maxc(qq'|roundqq');
    f = sumc(ln(chqq)) + captst*ln(lambada);
endif;

if kd==3;
     qq = -sumc(.5*ln(h0[ki:n,1]));
    qq = qq - sumc(.5*(nu+1)*ln(1+eta2[ki:n,1]./(h0[ki:n,1]*(nu-2))));
    f = captst*ln(lambada) + qq;
endif;

if kc==2;
  "log likelihood= ";; f;
  "Variance= ";;(z0[ki:n,1]'z0[ki:n,1])/(captst-1);
ess=(z0[ki:n,1]'z0[ki:n,1]);
  "Residuals Sum of Squares= ";;ess;
rsq=1-ess/((y[ki:n,1]-meanc(y[ki:n,1]))'*(y[ki:n,1]-meanc(y[ki:n,1])));
  "R-squared= ";;rsq;
aic=ln(ess/(captst-1))+ 2*(spar-1) /(captst-spar+1);
"AIC= ";;aic;
  "Mean of sqrt(h)= ";;meanc(sqrt(h0[ki:n,1]));
  "Mean of h= ";;meanc(h0[ki:n,1]);
  "";
  "Analisis of residuals";

yxx=z0[ki:n,1]./sqrt(h0[ki:n,1]);
"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;
"";
"Level testing:";
rho=zeros(kau,1);

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

k=1;
do until k > kau;
rho[k,1]=(yxx[1:captst-k,1]-meanc(yxx))'*(yxx[1+k:captst,1]-meanc(yxx))/captst;
rho[k,1]=rho[k,1]/vrho;
k=k+1;
endo;
"First kau autocorrelations:";
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);

"";
"Squared testing:";
yxx=(z0[ki:n,1]./sqrt(h0[ki:n,1]))^2;

rho=zeros(kau,1);

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

k=1;
do until k > kau;
rho[k,1]=(yxx[1:captst-k,1]-meanc(yxx))'*(yxx[1+k:captst,1]-meanc(yxx))/captst;
rho[k,1]=rho[k,1]/vrho;
k=k+1;
endo;
"First kau autocorrelations:";
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);

    output file=tango off;
    output file=tango.var reset;
/*
"OUTLIERS";
"  OBS      ERROR    VARIANCE    NON-LIN P-value";
*/
    it=ki;
    itl=0;
    ith=0;
    do until it >n;
    v0[it,1]=abs(z0[it,1])/sqrt(h0[it,1]);
if kvarp==1;
        if kd==1;
            prob[it,1]=cdfnc(v0[it,1]);
        endif;
        if kd==3;
            prob[it,1]=cdftc(v0[it,1],nu);
        endif;
if mnl ==1;
    y[it,1]~z0[it,1]~h0[it,1];;
 else;
    y[it,1]~z0[it,1]~h0[it,1];
endif;
        if mnl==1;
            mnf[it,1];
        endif;
endif;

        if nlarch==1;
            nf[it-ki+1,1];;
        endif;
/*
        if prob[it,1] > 0.975;
@            prob[it,1];;"  *";@
            prob[it,1];
            itl = itl+1;
        elseif prob[it,1] < 0.025;
            prob[it,1];
@            prob[it,1];;"  *"; @
            ith = ith+1;
        else;
            prob[it,1];
        endif;
*/
    it = it+1;
    endo;
        output file=tango.var off;

        output file=tango on;
/*
"";"";"Number of observations below .025 level:";itl;
      "Number of observations above .975 level:";ith;
*/

"";
"One step Variance forecast errors";
msx=z0[ki:n,1]^2-h0[ki:n,1];

"MSE:";;sumc(msx^2)/n;
"MAE:";;sumc(abs(msx))/n;
lsx=ln(z0[ki:n,1]^2)-ln(h0[ki:n,1]);
"ALES:";;sumc(lsx^2)/n;
"ALE:";;sumc(abs(lsx))/n;
"";
if garchb>0;
 "Persistence parameter:";;a1+b1+l1/2;
"Unconditional variance:";;a0/(1-(a1+b1+l1/2));
endif;
"Variance for last (n) observation:";;h0[n,1];
"";
if mk>0;
 z00 = zeros(mk,1);

  z00 = y[n+1:n+mk,1] - alpha0 ;
  if exog >0;
  z00 = z00 - xx[n+1:n+mk,.]*beta;
  endif;
    iz = 1;
    do until iz > pphi;
       z00 = z00 - phi[iz,1]*y[n:n+mk-1,1];
    iz = iz+1;
    endo;

"Nth Step Ahead Forecast - Squared Forecast Error - Absolute Forecast Error";
  j=1;
  sab=0;
  hf=zeros(mk,1);
  do  until j > mk;
    sab=sab+(a1+b1+l1/2)^(j-1);
    hf[j,1]=a0*sab + ((a1+b1+l1/2)^j)*h0[n,1];
    serr=(z00[j,1]^2-hf[j,1])^2;
    aerr=abs(z00[j,1]^2-hf[j,1]);
    j~hf[j,1]~serr~aerr;
    j=j+1;
  endo;
  endif;
rsq~aic;
endif;

retp(-f);
endp;
/* ----------------------------------------------------------- */