@ PROC GARCH @ @ This proc will estimate a GARCH(p,q) model @ @ for the variance and a user-specified model @ @ for the mean of a series. @ @ @ @ INPUT: Tx1 dependent variable @ @ TxK independent variables @ @ number of lags of resids^2 in GARCH model @ @ number of lags of conditional variance @ @ limits: 1=calculate Wald CL's 0 = do not @ @ startvalues = vector of inital values @ @ OUTPUT: @ @ b, a (K+p+q+1)x1 vector of parameters @ @ f, scalar, function at minimum @ @ g, (K+p+q+1)x1 vector, gradient at min @ @ c11, Wald confidence limits @ @ ret, scalar, CML return code @ @ hhat, Tx1 vector of conditional variance @ @ estimates at each point in time @ /* GARCH Code found at: http://faculty.washington.edu/rons/garch.html A site connected to Ronald Schoenberg's site Modified by Andrew Patton email: apatton@econ.ucsd.edu Fri, Sep 3, 1999. This code is provided gratis, without performance warrantee. */ library cml; cmlset; __output=0; _p=0;_q=0;_numx=0; @ setting these up for future reference @ proc (6) = garch(depvar,indepvar,garchp,garchq,limits,startvalues); local p,q,z,numx,r,start,b,f,g,vcv,ret,cl1,_ww_, y, _numx, _p, _q, u2, hhat, u; p = garchp; q = garchq; y = varput(p,"_p"); @ variable made global so it can be used in other proc's @ y = varput(q,"_q"); @ variable made global so it can be used in other proc's @ /* 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. */ if indepvar/=0; z = depvar~indepvar; numx = cols(indepvar); _cml_ParNames = ("indep"$+ftocv(seqa(1,1,cols(indepvar)),1,0)) | "kappa" | (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) | (0$+"alpha"$+ftocv(seqa(1,1,q),1,0)); else; z = depvar; numx = 0; _cml_ParNames = "omega" | (0$+"beta"$+ftocv(seqa(1,1,p),1,0)) | (0$+"alpha"$+ftocv(seqa(1,1,q),1,0)); endif; y = varput(numx,"_numx"); @ variable made global so it can be used in other proc's @ /* constraints ** kappa > 0.0001, delta .>= 0, alpha .>= 0, sumc(delta|alpha) < 1 */ _ww_ = { -1e250 1e250 }; _cml_Bounds = ones(numx+1+p+q,2).*_ww_; _cml_Bounds[numx+1,1] = .0001; _cml_Bounds[numx+2:numx+p+q+1,1] = zeros(p+q,1); _cml_c = zeros(1,rows(_cml_Parnames)); _cml_c[1,numx+2:numx+p+q+1] = -ones(1,p+q); _cml_d = -.9999; /* Constraints for IGARCH model, uncomment these lines and ** comment out the lines above for _cml_c and _cml_d ** ** _cml_a = zeros(1,rows(_cml_Parnames)); ** _cml_a[1,numx+2:numx+p+q+1] = ones(1,p+q); ** _cml_b = 1; */ /* Set weights for first maxc(p|q) observations to zero ** to keep them out of the log-likelihood */ r = maxc(p|q); __weight = (rows(z)/(rows(z)-r))*ones(rows(z),1); __weight[1:r] = zeros(r,1); @ start values @ if indepvar/=0; start = z[.,1] / z[.,2:cols(z)]; start = start | 1 | .01 * ones(p+q,1); else; start = startvalues; endif; _cml_GradProc = &grd; { b,f,g,vcv,ret } = cml( z, 0, &garchll, start ); if indepvar/=0; u = (depvar - indepvar * b[1:cols(indepvar)]); _numx = cols(indepvar); else; u = (depvar); _numx=0; endif; u2 = u.^2; hhat = recserar(b[_numx+1]+_shape(u2,garchq,1)*b[_numx+garchp+2:_numx+garchp+garchq+1], (meanc(u2)-meanc(u).^2)*ones(garchp,1),b[_numx+2:_numx+garchp+1]); @ initial value for hhat is unconditional error variance @ if limits==1; __title = "Wald confidence limits"; cl1 = cmlclimits(b,vcv); call cmlclprt(b,f,g,cl1,ret); else; cl1 = miss(0,0); @ print "parameter estimates" b; @ @ print "variance covariance matrix for parameter estimates" vcv; @ endif; retp(b,f,g,cl1,ret,hhat); endp; @ %%%%%%%%%% some supplementary procedures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @ proc garchll( b, x ); local u,u2,kappa,delta,alpha,h; if _numx>0; u = (x[.,1] - x[.,2:cols(x)] * b[1:_numx]); else; u = (x[.,1]); endif; u2 = u.^2; kappa = b[_numx+1]; delta = b[_numx+2:_numx+_p+1]; alpha = b[_numx+_p+2:_numx+_p+_q+1]; h = recserar(kappa+_shape(u2,_q,1)*alpha,(meanc(u2)-meanc(u).^2)*ones(_p,1),delta); retp( -0.5*( (u2 ./ h) + ln(2 * pi) + ln(h) ) ); endp; proc grd(b,x); @ original code for the gradient @ local j,u,u2,kappa,delta,alpha,r,h,s,g; local p,q,numx; p = _p; q = _q; numx = _numx; u = x[.,1] - x[.,2:cols(x)] * b[1:numx]; u2 = u .* u; kappa = b[numx+1]; delta = b[numx+2:numx+p+1]; alpha = b[numx+p+2:numx+p+q+1]; g = ones(rows(x),rows(b)); s = _shape(u2,q,1); h = recserar(kappa+s*alpha,(meanc(u2)-meanc(u).^2)*ones(p,1),delta); g[.,1:numx] = -2*_dshape(u.*x[.,2:cols(x)],q,alpha); g[.,numx+2:numx+p+1] = _shape(h,p,-1); g[.,numx+p+2:numx+p+q+1] = s; g = (0.5*(u2./h - 1)./h) .* recserar(g,zeros(p,cols(g)),delta.*ones(1,cols(g))); g[.,1:numx] = g[.,1:numx] + (u./h).*x[.,2:cols(x)]; retp(g); endp; @ Modified gradient code, trying to deal with the case where @ @ there is NOTHING in the model for the mean. Useful for @ @ estimating autoregressive conditional duration (ACD) models @ proc grd2(b,x); local j,u,u2,kappa,delta,alpha,r,h,s,g; if _numx>0; u = x[.,1] - x[.,2:cols(x)] * b[1:_numx]; else; u = x[.,1]; endif; u2 = u .* u; kappa = b[_numx+1]; delta = b[_numx+2:_numx+_p+1]; alpha = b[_numx+_p+2:_numx+_p+_q+1]; g = ones(rows(x),rows(b)); s = _shape(u2,_q,1); h = recserar(kappa+s*alpha,(meanc(u2)-meanc(u).^2)*ones(_p,1),delta); if _numx>0; g[.,1:_numx] = -2*_dshape(u.*x[.,2:cols(x)],_q,alpha); endif; g[.,_numx+2:_numx+_p+1] = _shape(h,_p,-1); g[.,_numx+_p+2:_numx+_p+_q+1] = s; g = (0.5*(u2./h - 1)./h) .* recserar(g,zeros(_p,cols(g)),delta.*ones(1,cols(g))); if _numx>0; g[.,1:_numx] = g[.,1:_numx] + (u./h).*x[.,2:cols(x)]; endif; retp(g); 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 _dshape(z,m,a); local y,j,n,v; n = rows(z) - m - 1; y = zeros(rows(z),cols(z)); v = seqa(1,1,m)' + seqa(0,1,n+1); j = 1; do until j > cols(z); y[m+1:rows(z),j] = reshape(z[v,j],n+1,m) * a; j = j + 1; endo; retp(y); endp;