/*
**  SGARCH model   (Schoenberg Dec95)
**
**  Estimates parameters of the SGARCH model.  An example and the data
**  set, as well as a discussion of the SGARCH model, can be found at
**  http://weber.u.washington.edu/~rons
**
** Note that the log-likelihood proc used in the first
** call to CML to estimate the mean and alpha for the
** first R observations can also be used in general for
** regression estimates based on the symmetric stable
** distribution.
**
**  This program was written by Ron Schoenberg, rons@u.washington.edu,
**  and is in the public domain.  While the author disclaims any
**  responsibility for the performance of this software, he would
**  appreciate receiving any comments.
**
**  This program requires GAUSS 3.2.7+ and the Constrained Maximum
**  Likelihood application module.
*/


  library cml;
  #include cml.ext;
  cmlset;


/* set order here for a GARCH(p,q) model */

  p = 1;
  q = 1;


/*
**  Choose number of observations to be used in
**  computing the initial scale of the error term.
**  These observations will not be included in the
**  SGARCH estimation
*/

   R = 20;


/* data  */

  dsn = "garch";

  depvar = { INDEX };
  indvar = { };


  output file = sgarch.out reset;

/*
**  Usually it will not be necessary to change any of the code
**  below this line.  Exceptions might be to specify novel
**  constraint configurations.
*/



/*
**  The dataset is read in using loadd because the GARCH log-likelihood
**  requires the entire data set and passing the data in as a matrix
**  will ensure that the complete data set will be passed to the
**  log-likelihood proc at each call to the function.
*/

  z = loadd(dsn);

  if not scalmiss(indvar);
     { depvar,dv,indvar,iv } = indices2(dsn,depvar,indvar);
  else;
    { depvar,dv } = indices(dsn,depvar);
  endif;

  if scalmiss(dv);
     errorlog "error: variable not found in dataset "$+dsn;
     end;
  endif;


  if not scalmiss(indvar);
    z = z[.,dv]~ones(rows(z),1)~z[.,iv];
    numx = rows(indvar) + 1;
  else;
    z = z[.,dv]~ones(rows(z),1);
    numx = 1;
  endif;



/*
** Globals for stable density procedure
*/

    clear s_,c5i_,a_, alf2_, alf1_, pd_, p_, pdd_, oldalf_;
    oldalf_ = -999.99;


/* this s is the gij of the writeup */

    s_ = { 1.8514190959e+2 -4.6769332663e+2 4.8424720302e+2 -1.7639153404e+2,
        -3.0236552164e+2 7.6351931975e+2 -7.8560342101e+2 2.8426313374e+2,
         4.4078923600e+2 -1.1181138121e+3 1.1548311335e+3 -4.1969666223e+2,
        -5.2448142165e+2 1.3224487717e+3 -1.3555648053e+3 4.8834079950e+2,
         5.3530435018e+2 -1.3374570340e+3 1.3660140118e+3 -4.9286099583e+2,
        -4.8988957866e+2 1.2091418165e+3 -1.2285872257e+3 4.4063174114e+2,
         3.2905528742e+2 -7.3211767697e+2 6.8183641829e+2 -2.2824291084e+2,
        -2.1495402244e+2 3.9694906604e+2 -3.3695710692e+2 1.0905855709e+2,
         2.1112581866e+2 -2.7921107017e+2 1.1717966020e+2 3.4394664342e+0,
        -2.6486798043e+2 1.1999093707e+2 2.1044841328e+2 -1.5110881541e+2,
         9.4105784123e+2 -1.7221988478e+3 1.4087544698e+3 -4.2472511892e+2,
        -2.1990475933e+3 4.2637720422e+3 -3.4723981786e+3 1.0174373627e+3,
         3.1047490290e+3 -5.4204210990e+3 4.2221052925e+3 -1.2345971177e+3,
        -5.1408260668e+3 1.1090264364e+4 -1.0270337246e+4 3.4243449595e+3,
         1.1215157876e+4 -2.4243529825e+4 2.1536057267e+4 -6.8490996103e+3,
        -1.8120631586e+4 3.1430132257e+4 -2.4164285641e+4 6.9126862826e+3,
         1.7388413126e+4 -2.2108397686e+4 1.3397999271e+4 -3.1246611987e+3,
        -7.2435775303e+3 4.3545399418e+3 2.3616155949e+2 -7.6571653073e+2,
        -8.7376725439e+3 1.5510852129e+4 -1.3789764138e+4 4.6387417712e+3 };

    c5i_ = { 1, 5, 10, 10, 5, 1 };


/*
** start values
*/

  if not scalmiss(indvar);
    start0 = z[.,1] / z[.,2:cols(z)];
  else;
    start0 = meanc(z[.,1]);
  endif;


/*
**  ESTIMATION OF INITIAL CONDITIONS
**
**  The first R observations are fit to a univariate
**  stable distribution.  The estimate of alpha here
**  will be used to initialize the error term sequence
*/

  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(numx+2,2).*_ww_;   /* set all to (-inf, +inf) */
  _cml_Bounds[numx+1,1] = .84;          /* alpha >= .84 */
  _cml_Bounds[numx+1,2] = 2;            /* alpha <= 2 */
  _cml_Bounds[numx+2,1] = .0001;        /* variance > .0001 */

  start = start0 | 1.8 | 1;


  if not scalmiss(indvar);
    _cml_ParNames = "CONST" | indvar | "alpha" | "h0";
  else;
    _cml_ParNames = "CONST" | "alpha" | "h0";
  endif;
  __output = 0;

  __title = "Estimate of alpha for initial R observations";
  { b,f,g,h,ret } = cml( z[1:r,.], 0, &stable, start );


  if ret /= 0;
      print "estimate of initial conditions with unknown";
      print "alpha failed - scale will now be estimated ";
      print "assuming alpha = 1.8";

      start = start0 | 1.8 | 1;
     _cml_Active = ones(numx+2,1);
     _cml_Active[numx+1] = 0;
     { b,f,g,h,ret } = cml( z[1:r,.], 0, &stable, start );
   endif;

   call cmlprt(b,f,g,h,ret);
    _h0_ = b[numx+2]^b[numx+1];

    cmlset;  @ reset globals @




/*
**  SGARCH ESTIMATION
**
**  constraints
**
**  kappa >= 0.001,
**  delta .>= 0,
**  theta .>= 0,
**  sumc(delta|theta) <= .999,
**  .84 <= alpha <= 2
**
**  alpha constrained to be greater than or equal to .84 because numerical
**  computation of stable density is not adequate for alpha < .84
**
**  alpha = 2 is Normal density, and alpha > 2 is undefined
*/

  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(numx+p+q+2,2).*_ww_;    /* set all to (-inf, +inf) */
  _cml_Bounds[numx+1,1] = .001;              /* kappa > 0 */
  _cml_Bounds[numx+2:numx+p+q+1,1] = zeros(p+q,1);  /* delta,theta >= 0 */
  _cml_Bounds[numx+p+q+2,1] = .84;           /* alpha >= .84 */
  _cml_Bounds[numx+p+q+2,2] = 2;             /* alpha <= 2 */

  _cml_c = zeros(1,rows(_cml_Bounds));       /* sumc(delta|theta) <= .999 */
  _cml_c[1,numx+2:numx+p+q+1] = -ones(1,p+q);
  _cml_d = -.999;


  start = start0 | 1 | .01 * ones(p+q,1) | 1.8;

   if not scalmiss(indvar);
    _cml_ParNames = "CONST" | indvar | "kappa" |
        (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) |
        (0$+"theta"$+ftocv(seqa(1,1,q),1,0)) | "alpha";
  else;
    _cml_ParNames = "CONST" | "kappa" |
        (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) |
        (0$+"theta"$+ftocv(seqa(1,1,q),1,0)) | "alpha";
  endif;



/*
**  Set weights for first R observations to zero
**  to keep them out of the log-likelihood
*/

  __weight = (rows(z)/(rows(z)-R))*ones(rows(z),1);
  __weight[1:R] = zeros(R,1);


/*
**  Call CML to generate maximum Likelihood estimates
*/

  __output = 20;

  { b,f,g,h,ret } = CML( z, 0, &sgarch, start );


/*
**  Statistical inference
*/

  if cond(h) < 1e16;

      __title = "Ordinary Wald confidence limits";
      clt = cmltlimits(b,h);
      call cmlclprt(b,f,g,clt,ret);

      __title = "Inequality constrained Wald confidence limits";
      clw = cmlclimits(b,h);
      call cmlclprt(b,f,g,clw,ret);
  else;

      print "covariance matrix of parameters not sufficiently"\
             " well conditioned to compute Wald confidence limits";
  endif;

  _cml_BootFname = "sboot";
  _cml_NumSample = 500;

  __title = "bootstrap confidence limits";
  { b,f,g,h,ret } = cmlboot( z, 0, &sgarch, start );

  clb = cmlblimits("sboot");
  call cmlclprt(b,f,g,clb,ret);

  output off;


/*************** procedures *******************/

proc stable( b, x );
  local u,h,alpha;

  u = x[.,1] - x[.,2:cols(x)] * b[1:numx];
  alpha = b[numx+1];
  h = b[numx+2];

  retp(ln(symstb1(u./h^(1/alpha),alpha))-ln(h)/alpha);

endp;



proc sgarch( b, x );
  local u,u2,kappa,delta,theta,alpha,h;

  kappa = b[numx+1];
  delta = b[numx+2:numx+p+1];
  theta = b[numx+p+2:numx+p+q+1];
  alpha = b[numx+p+q+2];

  u = x[.,1] - x[.,2:cols(x)] * b[1:numx];

  h = recserar(kappa+_shape(abs(u)^alpha,q,1)*theta,_h0_*ones(p,1),delta);

  retp(ln(symstb1(u./h^(1/alpha),alpha))-ln(h)/alpha);

endp;


proc _shape(z,m,r);
   local y,n,v;
   n = rows(z) - m - 1;

   y = zeros(rows(z),m);
   if r == 1;
       v = seqa(1,1,m)' + seqa(0,1,n+1);
   else;
       v = seqa(m,-1,m)' + seqa(0,1,n+1);
   endif;
   y[m+1:rows(z),.] = reshape(z[v],n+1,m);

   retp(y);
endp;



proc symstb1(x,alpha);

/*
** June 1993: GAUSS version
** If you use this program, please cite J. Huston McCulloch, "Numerical
** Approximation of the Symmetric Stable Distribution and Density," OSU
** Econ Dept., Oct. 1994.  x may be a column vector; alpha is a scalar
** returns sden = symmetric stable density
*/

    local r, q, znot, zn4, zn5;
    local i, i1, j, j1;
    local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
    local rp0, rpp0, rppp0, rp1, b, c;
    local xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden;
    local x11, x2a1, x2pp, cd, gd;

    /*
    if alpha > 2.000001;
        retp(error(0));
    endif;
    */

    if alpha >= 2;
        retp(exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) /
                                              3.5449077018110318);
    endif;

    if not (alpha $== oldalf_);

        znot = seqa(1,1,19) /20;
        zn4 = (1-znot)^4;
        a_ = 2^(1/alpha) - 1;
        alf2_ = 2 - alpha;
        alf1_ = alpha - 1;
        pialf = pi * alpha;
        sp0 = gamma(1/alpha) / pialf;
        xp0 = 1/(alpha * a_);
        xpp0 = xp0 * (1+alpha) / alpha;
        q = zeros(6,1);
        q[2] = -sp0 * xp0 + alf2_ * .318309881837907 + alf1_ *
                   .3405185360876554;
        q[3] = .5*(-sp0 * xpp0 + alf2_ * .6366197723675814 +
                            alf1_ * .5107778041314830);
        q[4] = (-sp0 * xpp0 * (1+2*alpha)/alpha + (gamma(3/alpha) / pialf) *
            xp0^3 + alf2_ * 1.2732395447351628 + alf1_ * 1.0288585763021882) /
 6;
        r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
        b = -sumc(q[2:4]) - sumc(r.*zn4.*(1-znot));
        c = -(a_^alpha) * gamma(alpha) * sin(pialf / 2) / pi + alf2_ / pi
                      - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
        q[5] = 5*b-c;
        q[6] = c-4*b;
        p_ = zeros(6,20);
        i = 0;
        do until i > 5;
            i1 = i + 1;
            j = 1;
            do until j > 19;
                j1 = j + 1;
                p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
                j = j1;
            endo;
            i = i1;
        endo;
        p_ = q + c5i_ .* p_;
        pd_ = p_[2:6,.] .* seqa (1,1,5);
        oldalf_ = alpha;
    endif;

    xa1 = 1 + a_ * abs(x);
    xa1a = xa1 ^ (-alpha);
    z = 1 - xa1a;
    zp = alpha * a_ * xa1a ./ xa1;
    x11 = 1/xa1a;
    x1 = x11 - 1;
    x1p = x11 .* x11;
    x2a1 = 1/sqrt(xa1a);
    x2 = (x2a1 -1) /.4142135623730951;
    x2p = x2a1 ./ (.8284271247461902 * xa1a);

    cd = 1 ./ (pi*(1 + x1.*x1));
    gd = exp(maxc(-700*ones(1,rows(x2)) | -(x2.*x2)'/4)) /
                                             3.5449077018110318;

    j1 = trunc(20 * z) + 1;
    rpz = (((pd_[5,j1]'.*z + pd_[4,j1]').*z + pd_[3,j1]').*z +
                              pd_[2,j1]').*z + pd_[1,j1]';
    sden = (alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz) .* zp;

    retp (sden);
endp;
