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