dlibrary recurs; library cml; #include cml.ext; cmlset; dsn = "garch"; depvar = { INDEX }; indvar = { }; /* set order here for a GARCH-M(p) model */ p = 1; output file = garchm.out reset; /* ** 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. */ 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; _cml_ParNames = "CONST" | indvar | "kappa" | (0$+"theta"$+ftocv(seqa(1,1,p),1,0)) | "gamma"; else; z = z[.,dv]~ones(rows(z),1); numx = 1; _cml_ParNames = "CONST" | "kappa" | (0$+"theta"$+ftocv(seqa(1,1,p),1,0)) | "gamma"; endif; /* ** constraints ** ** kappa > 0.01, theta .>= 0, sumc(theta) < 1 ** */ _ww_ = { -100 100 }; _cml_Bounds = ones(numx+2+p,2).*_ww_; _cml_Bounds[numx+1,1] = .01; _cml_Bounds[numx+2:numx+p+1,1] = zeros(p,1); _cml_c = zeros(1,rows(_cml_Parnames)); _cml_c[1,numx+2:numx+p+1] = -ones(1,p); _cml_d = -.99; /* ** Set weights for first p observations to zero ** to keep them out of the log-likelihood */ __weight = (rows(z)/(rows(z)-p))*ones(rows(z),1); __weight[1:p] = zeros(p,1); /* ** start values */ if not scalmiss(indvar); start = z[.,1] / z[.,2:cols(z)]; else; start = meanc(z[.,1]); endif; start = start | 1.1 | .01 * ones(p,1) | .001; _cml_Options = { XPROD }; { b,f,g,h,ret } = cml( z, 0, &garchm, start ); __title = "Wald confidence limits"; cl1 = cmlclimits(b,h); call cmlclprt(b,f,g,cl1,ret); __title = "LR confidence limits"; cl1 = cmlpflclimits(b,f,z,0,&garchm); call cmlclprt(b,f,g,cl1,ret); _cml_BootFname = "gmboot"; _cml_NumSample = 500; { b,f,g,h,ret } = cmlboot( z, 0, &garchm, start ); __title = "bootstrap confidence limits"; clb = cmlblimits("gmboot"); call cmlclprt(b,f,g,clb,ret); output off; /*************** procedures *******************/ proc garchm( b, x ); local u,beta,gamm,kappa,delta,theta,h,i,rx,cx; beta = b[1:numx]; kappa = b[numx+1]; theta = b[numx+2:numx+p+1]; gamm = b[numx+p+2]; rx = rows(x); cx = cols(x); u = zeros(rx,1); h = zeros(rx,1); h[1:p] = meanc((x[.,1] - x[.,2:cols(x)] * beta)^2) * ones(p,1); u[1:p] = h[1:p]; dllcall recursdll(h,u,x,rx,cx,p,kappa,theta,beta,gamm); /* ** recursdll does the following calculation: ** ** i = p + 1; ** do until i > rx; ** h[i] = kappa + h[i-p:i-1]'*theta; ** u[i] = x[i,1] - x[i,2:cx]' * beta - gamm * h[i]; ** i = i + 1; ** endo; */ retp( -0.5*( (u .* u) ./ h + ln(2 * pi) + ln(h) ) ); endp; /* ** Source code for recursdll: ** ** To compile use the following statement: ** ** cc -G librecurs.c -o librecurs.so -lm ** ** The copy the resulting library, i.e., librecurs.so, to ** gauss/dlib. ** ** ** int recursdll( double *h, double *u, double * x, double *rx, ** double *cx, double *p, double *kappa, double *theta, ** double *beta, double *gamm ) ** { ** int i, j, k, p0, rx0, cx0; ** double e; ** ** p0 = *p; ** rx0 = *rx; ** cx0 = *cx; ** ** for (i = p0; i < rx0; i++) { ** ** h[i] = *kappa; ** ** for (j = 0, k = i - 1; j < p0; j++, k--) { ** h[i] += theta[j] * h[k]; ** } ** ** u[i] = x[i*cx0] - (*gamm) * h[i]; ** for (j = 0, k = i*cx0 + 1; j < cx0 - 1; j++, k+=cx0) { ** u[i] -= beta[j] * x[k]; ** } ** } ** return 0; ** } */