1 /* Copyright (C) 1997-1999 Adrian Trapletti
2
3 This library is free software; you can redistribute it and/or
4 modify it under the terms of the GNU Library General Public
5 License as published by the Free Software Foundation; either
6 version 2 of the License, or (at your option) any later version.
7
8 This library is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 Library General Public License for more details.
12
13 You should have received a copy of the GNU Library General Public
14 License along with this library; if not, write to the Free
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 GARCH estimation
18
19 Reference: T. Bollerslev (1986): Generalized Autoregressive Conditional
20 Heteroscedasticity, Journal of Econometrics 31, 307-327. */
21
22
23 #include <R.h>
24
25
26 extern void F77_NAME(dsumsl) ();
27 extern void F77_NAME(dsmsno) ();
28 extern void F77_NAME(ddeflt) ();
29
30
31 #define BIG 1.0e+10 /* function value if the parameters are invalid */
32
33
34 static double dsqrarg;
35 #define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
36
37 static double dmaxarg1,dmaxarg2;
38 #define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
39 (dmaxarg1) : (dmaxarg2))
40
41
42 struct garch_handler /* used to set up the additional parameters used in calcf and calcg */
43 {
44 double* y; /* the time series to fit */
45 double* h; /* the conditional variance (cv) */
46 double* dh; /* dh_i/dp_j */
47 int n; /* the number of observations */
48 int p, q; /* GARCH(p,q) */
49 };
50
51 static struct garch_handler garch_h;
52
53
F77_SUB(calcf)54 static void F77_SUB(calcf) (int *pq, double *p, int *nf, double *f, int *uiparm,
55 double *urparm, void (*F77_SUB(ufparm))(void))
56 /* compute negative log likelihood apart from the constant and the pre-sample values */
57 {
58 int i, j, ok;
59 int maxpq = DMAX(garch_h.p,garch_h.q);
60 double temp = 0.0;
61 double sum = 0.0;
62
63 ok = 1;
64 if (p[0] <= 0.0) ok = 0;
65 for (i=1; i<(*pq); i++)
66 if (p[i] < 0.0) ok = 0;
67 if (ok) /* parameters are valid */
68 {
69 for (i=maxpq; i<garch_h.n; i++) /* loop over time */
70 { /* compute cv at time i */
71 temp = p[0]; /* compute ARCH part of cv */
72 for (j=1; j<=garch_h.q; j++)
73 temp += p[j]*DSQR(garch_h.y[i-j]);
74 for (j=1; j<=garch_h.p; j++) /* compute GARCH part of cv */
75 temp += p[garch_h.q+j]*garch_h.h[i-j];
76 sum += log(temp)+DSQR(garch_h.y[i])/temp; /* compute eq. 18 */
77 garch_h.h[i] = temp; /* assign cv at time i */
78 }
79 (*f) = 0.5*sum;
80 }
81 else /* parameters are invalid */
82 (*f) = BIG;
83 }
84
F77_SUB(calcg)85 static void F77_SUB(calcg) (int *pq, double *p, int *nf, double *dp, int *uiparm,
86 double *urparm, void (*F77_SUB(ufparm))(void))
87 /* compute derivative of negative log likelihood */
88 {
89 int i, j, k;
90 int maxpq = DMAX(garch_h.p,garch_h.q);
91 double temp1, temp2, temp3;
92
93 for (k=0; k<(*pq); k++) /* initialize */
94 dp[k] = 0.0;
95 for (i=maxpq; i<garch_h.n; i++) /* loop over time */
96 { /* compute cv at time i and derivatives dh_i/dp_j */
97 temp1 = p[0]; /* compute ARCH part of cv */
98 for (j=1; j<=garch_h.q; j++)
99 temp1 += p[j]*DSQR(garch_h.y[i-j]);
100 for (j=1; j<=garch_h.p; j++) /* compute GARCH part of cv */
101 temp1 += p[garch_h.q+j]*garch_h.h[i-j];
102 garch_h.h[i] = temp1; /* assign cv at time i */
103 temp2 = 0.5*(1.0-DSQR(garch_h.y[i])/temp1)/temp1; /* compute dl_i/dh_i, eq. 19 */
104 temp3 = 1.0; /* compute dh_i/dp_0, eq. 21 */
105 for (j=1; j<=garch_h.p; j++)
106 temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)];
107 garch_h.dh[(*pq)*i] = temp3; /* assign dh_i/dp_0 */
108 dp[0] += temp2*temp3; /* assign dl_i/dp_0 = dl_i/dh_i * dh_i/dp_0 */
109 for (k=1; k<=garch_h.q; k++) /* compute dl_i/dp_k for the ARCH part, eq. 19 */
110 {
111 temp3 = DSQR(garch_h.y[i-k]); /* compute dh_i/dp_k, eq. 21 */
112 for (j=1; j<=garch_h.p; j++)
113 temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)+k];
114 garch_h.dh[(*pq)*i+k] = temp3; /* assign dh_i/dp_k */
115 dp[k] += temp2*temp3; /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
116 }
117 for (k=1; k<=garch_h.p; k++) /* compute dl_i/dp_k for the GARCH part, eq. 19 */
118 {
119 temp3 = garch_h.h[i-k]; /* compute dh_i/dp_k, eq. 21 */
120 for (j=1; j<=garch_h.p; j++)
121 temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)+garch_h.q+k];
122 garch_h.dh[(*pq)*i+garch_h.q+k] = temp3; /* assign dh_i/dp_k */
123 dp[garch_h.q+k] += temp2*temp3; /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
124 }
125 }
126 }
127
F77_SUB(ufparm)128 static void F77_SUB(ufparm) ()
129 {
130 error ("fatal error in fit_garch ()\n");
131 }
132
133
tseries_fit_garch(double * y,int * n,double * par,int * p,int * q,int * itmax,double * afctol,double * rfctol,double * xctol,double * xftol,double * fret,int * agrad,int * trace)134 void tseries_fit_garch (double *y, int *n, double *par, int *p, int *q,
135 int *itmax, double *afctol, double *rfctol,
136 double *xctol, double *xftol, double *fret,
137 int *agrad, int *trace)
138 /* fit a GARCH (p, q) model
139
140 Input:
141
142 y[0..n-1] series to fit
143 par[0..p+q] initial parameter estimates
144 p, q model orders
145 itmax maximum number of iterations
146 afctol absolute function convergence tolerance
147 rfctol relative function convergence tolerance
148 xctol x-convergence tolerance
149 xftol false convergence tolerance
150 agrad estimation with analytical/numerical gradient
151 trace output yes/no
152
153 Output:
154
155 par[0..p+q] parameter estimates at minimum
156 fret function value at minimum
157 */
158 {
159 int i, j, pq, liv, lv, alg;
160 double *d, *v;
161 int *iv;
162 double var;
163
164 /* set up general optimizer parameters to default values */
165 pq = (*p)+(*q)+1;
166 d = Calloc (pq, double);
167 for (i=0; i<pq; i++)
168 d[i] = 1.0;
169 liv = 60;
170 iv = Calloc (liv, int);
171 lv = 77+pq*(pq+17)/2;
172 v = Calloc (lv, double);
173 alg = 2;
174 F77_CALL(ddeflt) (&alg, iv, &liv, &lv, v);
175 iv[0] = 12;
176
177 /* set up user defined optimizer parameters */
178 iv[16] = 2*(*itmax);
179 iv[17] = (*itmax);
180 if (*trace)
181 iv[20] = 6;
182 else
183 iv[20] = 0;
184 v[30] = (*afctol);
185 v[31] = (*rfctol);
186 v[32] = (*xctol);
187 v[33] = (*xftol);
188
189 /* set handler values */
190 garch_h.p = (*p); garch_h.q = (*q); garch_h.n = (*n); garch_h.y = y;
191 garch_h.h = Calloc ((*n), double);
192 garch_h.dh = Calloc ((*n)*pq, double);
193 var = 0.0;
194 for (i=0; i<(*n); i++) /* estimate unconditional variance (uv) */
195 var += DSQR(y[i]);
196 var /= (double) (*n);
197 for (i=0; i<DMAX((*p),(*q)); i++) /* initialize */
198 {
199 garch_h.h[i] = var; /* with uv */
200 garch_h.dh[pq*i] = 1.0; /* dh_i/dp_0 with 1 */
201 for (j=1; j<pq; j++)
202 garch_h.dh[pq*i+j] = 0.0; /* dh_i/dp_j with 0 */
203 }
204
205 if (*agrad) /* estimation with analytical gradient */
206 {
207 if (*trace) Rprintf ("\n ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** \n\n");
208 F77_CALL(dsumsl) (&pq, d, par, F77_SUB(calcf), F77_SUB(calcg),
209 iv, &liv, &lv, v, NULL, NULL, F77_SUB(ufparm));
210 if (*trace) Rprintf ("\n");
211 }
212 else /* estimation with numerical gradient */
213 {
214 if (*trace) Rprintf ("\n ***** ESTIMATION WITH NUMERICAL GRADIENT ***** \n\n");
215 F77_CALL(dsmsno) (&pq, d, par, F77_SUB(calcf), iv,
216 &liv, &lv, v, NULL, NULL, F77_SUB(ufparm));
217 if (*trace) Rprintf ("\n");
218 }
219
220 /* return function value */
221 (*fret) = v[9];
222
223 /* free memory */
224 Free (d);
225 Free (iv);
226 Free (v);
227 Free (garch_h.h);
228 Free (garch_h.dh);
229 }
230
231
tseries_pred_garch(double * y,double * h,int * n,double * par,int * p,int * q,int * genuine)232 void tseries_pred_garch (double *y, double *h, int *n, double *par,
233 int *p, int *q, int *genuine)
234 /* predict cv with a GARCH (p, q) model
235
236 Input:
237
238 y[0..n-1] series to predict
239 par[0..p+q] parameters of the GARCH (p, q)
240 p, q model orders
241 genuine logical indicating if a genuine prediction is computed
242
243 Output:
244
245 h[0..N] predicted cv, where N = n for genuine prediction, and N = n-1 otherwise
246 */
247 {
248 double var, temp;
249 int i, j, maxpq, N;
250
251 if (*genuine) N = (*n)+1;
252 else N = (*n);
253 maxpq = DMAX((*p),(*q));
254 var = 0.0;
255 for (i=1; i<=(*p)+(*q); i++) /* compute uv */
256 var += par[i];
257 var = par[0]/(1.0-var);
258 for (i=0; i<maxpq; i++) /* initialize with uv */
259 h[i] = var;
260 for (i=maxpq; i<N; i++) /* loop over time */
261 { /* compute cv at time i */
262 temp = par[0]; /* compute ARCH part of cv */
263 for (j=1; j<=(*q); j++)
264 temp += par[j]*DSQR(y[i-j]);
265 for (j=1; j<=(*p); j++) /* compute GARCH part of cv */
266 temp += par[(*q)+j]*h[i-j];
267 h[i] = temp; /* assign cv at time i */
268 }
269 }
270
271
tseries_ophess_garch(double * y,int * n,double * par,double * he,int * p,int * q)272 void tseries_ophess_garch (double *y, int *n, double *par, double *he,
273 int *p, int *q)
274 /* Compute outer product approximation of the hessian of the
275 negative log likelihood of a GARCH (p, q) model at given parameter
276 estimates
277
278 Input:
279
280 y[0..n-1] time series
281 par[0..p+q] parameter estimates of the GARCH (p, q)
282 p, q model orders
283
284 Output:
285
286 he[0..(p+q+1)*(p+q+1)-1] predicted cv
287 */
288 {
289 double var, temp1, temp2, temp3;
290 int i, j, k, pq;
291 double *h, *dh, *dpar;
292
293 pq = (*p)+(*q)+1;
294 h = Calloc ((*n), double);
295 dh = Calloc ((*n)*pq, double);
296 dpar = Calloc (pq, double);
297 var = 0.0;
298 for (i=0; i<(*n); i++) /* estimate uv */
299 var += DSQR(y[i]);
300 var /= (double) (*n);
301 for (i=0; i<DMAX((*p),(*q)); i++) /* initialize */
302 {
303 h[i] = var; /* with uv */
304 dh[pq*i] = 1.0; /* dh_i/dp_0 with 1 */
305 for (j=1; j<pq; j++)
306 dh[pq*i+j] = 0.0; /* dh_i/dp_j with 0 */
307 }
308 for (k=0; k<pq; k++) /* initialize */
309 for (j=0; j<pq; j++)
310 he[pq*k+j] = 0.0;
311 for (i=DMAX((*p),(*q)); i<(*n); i++) /* loop over time */
312 { /* compute cv at time i and derivatives dh_i/dp_j */
313 temp1 = par[0]; /* compute ARCH part of cv */
314 for (j=1; j<=(*q); j++)
315 temp1 += par[j]*DSQR(y[i-j]);
316 for (j=1; j<=(*p); j++) /* compute GARCH part of cv */
317 temp1 += par[(*q)+j]*h[i-j];
318 h[i] = temp1; /* assign cv at time i */
319 temp2 = 0.5*(1.0-DSQR(y[i])/temp1)/temp1; /* compute dl_i/dh_i, eq. 19 */
320 temp3 = 1.0; /* compute dh_i/dp_0, eq. 21 */
321 for (j=1; j<=(*p); j++)
322 temp3 += par[(*q)+j]*dh[pq*(i-j)];
323 dh[pq*i] = temp3; /* assign dh_i/dp_0 */
324 dpar[0] = temp2*temp3; /* assign dl_i/dp_0 = dl_i/dh_i * dh_i/dp_0 */
325 for (k=1; k<=(*q); k++) /* compute dl_i/dp_k for the ARCH part, eq. 19 */
326 {
327 temp3 = DSQR(y[i-k]); /* compute dh_i/dp_k, eq. 21 */
328 for (j=1; j<=(*p); j++)
329 temp3 += par[(*q)+j]*dh[pq*(i-j)+k];
330 dh[pq*i+k] = temp3; /* assign dh_i/dp_k */
331 dpar[k] = temp2*temp3; /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
332 }
333 for (k=1; k<=(*p); k++) /* compute dl_i/dp_k for the GARCH part, eq. 19 */
334 {
335 temp3 = h[i-k]; /* compute dh_i/dp_k, eq. 21 */
336 for (j=1; j<=(*p); j++)
337 temp3 += par[(*q)+j]*dh[pq*(i-j)+(*q)+k];
338 dh[pq*i+(*q)+k] = temp3; /* assign dh_i/dp_k */
339 dpar[(*q)+k] = temp2*temp3; /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
340 }
341 for (k=0; k<pq; k++) /* compute outer product approximation, p. 317 */
342 for (j=0; j<pq; j++)
343 he[pq*k+j] += dpar[k]*dpar[j];
344 }
345 Free (h);
346 Free (dh);
347 Free (dpar);
348 }
349
350