1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "version.h"
22
23 #include <glib.h>
24
25 /* Levin, Lin and Chu (Journal of Econometrics, 2002), Table 2:
26 correction factors for mean (mu) and standard deviation (s)
27 in context of panel unit-root statistic.
28
29 T = sample size (actually \tilde{T})
30
31 k = 1: no constant
32 k = 2: constant included
33 k = 3: constant plus trend
34 */
35
get_LLC_corrections(int T,int k,double * mu,double * sigma)36 static int get_LLC_corrections (int T, int k, double *mu, double *sigma)
37 {
38 const double LLCfac[] = {
39 /* T mu1 s1 mu2 s2 mu3 s3 */
40 25, 0.004, 1.049, -0.554, 0.919, -0.703, 1.003,
41 30, 0.003, 1.035, -0.546, 0.889, -0.674, 0.949,
42 35, 0.002, 1.027, -0.541, 0.867, -0.653, 0.906,
43 40, 0.002, 1.021, -0.537, 0.850, -0.637, 0.871,
44 45, 0.001, 1.017, -0.533, 0.837, -0.624, 0.842,
45 50, 0.001, 1.014, -0.531, 0.826, -0.614, 0.818,
46 60, 0.001, 1.011, -0.527, 0.810, -0.598, 0.780,
47 70, 0.000, 1.008, -0.524, 0.798, -0.587, 0.751,
48 80, 0.000, 1.007, -0.521, 0.789, -0.578, 0.728,
49 90, 0.000, 1.006, -0.520, 0.782, -0.571, 0.710,
50 100, 0.000, 1.005, -0.518, 0.776, -0.566, 0.695,
51 250, 0.000, 1.001, -0.509, 0.742, -0.533, 0.603,
52 0, 0.000, 1.000, -0.500, 0.707, -0.500, 0.500 /* \infty */
53 };
54 int c = 0, err = 0;
55
56 if (k > 0 && k < 4) {
57 c = 2 * k - 1;
58 } else {
59 err = E_DATA;
60 }
61
62 if (!err) {
63 int i, r = 12;
64
65 for (i=0; i<12; i++) {
66 if (T <= LLCfac[7*i]) {
67 r = i;
68 break;
69 }
70 }
71 *mu = LLCfac[7*r+c];
72 *sigma = LLCfac[7*r+c+1];
73 }
74
75 return err;
76 }
77
78 /* detrend \delta y for Levin-Lin-Chu case 3 */
79
LLC_detrend(gretl_matrix * dy)80 static int LLC_detrend (gretl_matrix *dy)
81 {
82 gretl_matrix *X, *b;
83 int t, T = dy->rows;
84 int err;
85
86 X = gretl_matrix_alloc(T, 2);
87 b = gretl_matrix_alloc(2, 1);
88
89 if (X == NULL || b == NULL) {
90 err = E_ALLOC;
91 } else {
92 for (t=0; t<T; t++) {
93 gretl_matrix_set(X, t, 0, 1.0);
94 gretl_matrix_set(X, t, 1, t+1);
95 }
96 err = gretl_matrix_ols(dy, X, b, NULL, NULL, NULL);
97 }
98
99 if (!err) {
100 for (t=0; t<T; t++) {
101 /* replace with detrended values */
102 dy->val[t] -= (b->val[0] + b->val[1] * (t+1));
103 }
104 }
105
106 gretl_matrix_free(X);
107 gretl_matrix_free(b);
108
109 return err;
110 }
111
112 /* long-run variance calculation as per LLC */
113
LLC_lrvar(gretl_matrix * vdy,int K,int m,int * err)114 static double LLC_lrvar (gretl_matrix *vdy, int K, int m, int *err)
115 {
116 double w, s21 = 0, s22 = 0;
117 double *dy = vdy->val;
118 int T = vdy->rows;
119 int t, j;
120
121 if (m == 3) {
122 /* subtract linear trend */
123 *err = LLC_detrend(vdy);
124 if (*err) {
125 return NADBL;
126 }
127 } else if (m == 2) {
128 /* subtract the mean */
129 double dybar = 0;
130
131 for (t=0; t<T; t++) {
132 dybar += dy[t];
133 }
134 dybar /= T;
135 for (t=0; t<T; t++) {
136 dy[t] -= dybar;
137 }
138 }
139
140 for (t=0; t<T; t++) {
141 s21 += dy[t] * dy[t];
142 }
143
144 for (j=1; j<=K; j++) {
145 w = 1.0 - j /((double) K + 1);
146 for (t=j; t<T; t++) {
147 s22 += w * dy[t] * dy[t-j];
148 }
149 }
150
151 return (s21 + 2 * s22) / T;
152 }
153
154 /* In case we got a list of individual-specific ADF order terms,
155 check that it makes sense and do some basic accounting.
156 */
157
LLC_check_plist(const int * list,int N,int * pmax,int * pmin,double * pbar)158 static int LLC_check_plist (const int *list, int N, int *pmax, int *pmin,
159 double *pbar)
160 {
161 int err = 0;
162
163 if (list == NULL || list[0] == 0) {
164 err = E_DATA;
165 } else if (list[0] > 1 && list[0] != N) {
166 err = E_DATA;
167 } else {
168 int i;
169
170 *pmax = *pmin = *pbar = 0;
171
172 for (i=1; i<=list[0]; i++) {
173 if (list[i] < 0) {
174 err = E_DATA;
175 break;
176 }
177 if (list[i] > *pmax) {
178 *pmax = list[i];
179 }
180 if (i == 1 || list[i] < *pmin) {
181 *pmin = list[i];
182 }
183 *pbar += list[i];
184 }
185
186 if (list[0] > 1) {
187 *pbar /= N;
188 }
189 }
190
191 return err;
192 }
193
LLC_sample_check(int N,int t1,int t2,int m,const int * plist,int * NT)194 static int LLC_sample_check (int N, int t1, int t2, int m,
195 const int *plist, int *NT)
196 {
197 int i, p, minT, Ti;
198 int err = 0;
199
200 *NT = 0;
201
202 for (i=1; i<=plist[0] && !err; i++) {
203 p = plist[i];
204 minT = m + p + 1; /* ensure df > 0 */
205 if (minT < 4) {
206 minT = 4;
207 }
208 /* Ti is the regression-usable series length, after
209 accounting for required lags
210 */
211 Ti = t2 - t1 + 1 - (1 + p);
212 if (Ti < minT) {
213 err = E_DATA;
214 } else if (plist[0] == 1) {
215 *NT = N * Ti;
216 } else {
217 *NT += Ti;
218 }
219 }
220
221 return err;
222 }
223
DF_test_spec(int m)224 static const char *DF_test_spec (int m)
225 {
226 const char *tests[] = {
227 N_("test without constant"),
228 N_("test with constant"),
229 N_("with constant and trend"),
230 };
231
232 if (m > 0 && m < 4) {
233 return tests[m-1];
234 } else {
235 return "";
236 }
237 }
238
239 #define LLC_DEBUG 0
240
241 /* Levin-Lin-Chu panel unit-root test */
242
real_levin_lin(int vnum,const int * plist,DATASET * dset,gretlopt opt,PRN * prn)243 int real_levin_lin (int vnum, const int *plist, DATASET *dset,
244 gretlopt opt, PRN *prn)
245 {
246 int verbose = opt & OPT_V;
247 int u0 = dset->t1 / dset->pd;
248 int uN = dset->t2 / dset->pd;
249 int N = uN - u0 + 1; /* units in sample range */
250 gretl_matrix_block *B;
251 gretl_matrix *y, *b;
252 gretl_matrix *dy, *X, *ui;
253 gretl_matrix *e, *ei, *v, *vi;
254 gretl_matrix *eps;
255 double pbar, SN = 0;
256 int t, t1, t2, T, NT;
257 int s, pt1, pt2, dyT;
258 int i, j, k, K, m;
259 int p, pmax, pmin;
260 int bigrow, p_varies = 0;
261 int err;
262
263 err = LLC_check_plist(plist, N, &pmax, &pmin, &pbar);
264
265 if (err) {
266 return err;
267 }
268
269 /* the 'case' (1 = no const, 2 = const, 3 = const + trend */
270 m = 2; /* the default */
271 if (opt & OPT_N) {
272 /* --nc */
273 m = 1;
274 } else if (opt & OPT_T) {
275 /* --ct */
276 m = 3;
277 }
278
279 /* does p vary by individual? */
280 if (pmax > pmin) {
281 p_varies = 1;
282 }
283 p = pmax;
284
285 /* the max number of regressors */
286 k = m + pmax;
287
288 t1 = t2 = 0;
289
290 /* check that we have a useable common sample */
291
292 for (i=0; i<N && !err; i++) {
293 int pt1 = (i + u0) * dset->pd;
294 int t1i, t2i;
295
296 dset->t1 = pt1;
297 dset->t2 = dset->t1 + dset->pd - 1;
298 err = series_adjust_sample(dset->Z[vnum], &dset->t1, &dset->t2);
299 t1i = dset->t1 - pt1;
300 t2i = dset->t2 - pt1;
301 if (i == 0) {
302 t1 = t1i;
303 t2 = t2i;
304 } else if (t1i != t1 || t2i != t2) {
305 err = E_MISSDATA;
306 break;
307 }
308 }
309
310 if (!err) {
311 err = LLC_sample_check(N, t1, t2, m, plist, &NT);
312 }
313
314 if (!err) {
315 int Tbar = NT / N;
316
317 /* the biggest T we'll need for regressions */
318 T = t2 - t1 + 1 - (1 + pmin);
319
320 /* Bartlett lag truncation (LLC, p. 14) */
321 K = (int) floor(3.21 * pow(Tbar, 1.0/3));
322 if (K > Tbar - 3) {
323 K = Tbar - 3;
324 }
325
326 /* full length of dy vector */
327 dyT = t2 - t1;
328
329 B = gretl_matrix_block_new(&y, T, 1,
330 &dy, dyT, 1,
331 &X, T, k,
332 &b, k, 1,
333 &ui, T, 1,
334 &ei, T, 1,
335 &vi, T, 1,
336 &e, NT, 1,
337 &v, NT, 1,
338 &eps, NT, 1,
339 NULL);
340 if (B == NULL) {
341 err = E_ALLOC;
342 }
343 }
344
345 if (err) {
346 return err;
347 }
348
349 if (m > 1) {
350 /* constant in first column, if wanted */
351 for (t=0; t<T; t++) {
352 gretl_matrix_set(X, t, 0, 1.0);
353 }
354 }
355
356 if (m == 3) {
357 /* trend in second column, if wanted */
358 for (t=0; t<T; t++) {
359 gretl_matrix_set(X, t, 1, t+1);
360 }
361 }
362
363 if (verbose) {
364 pputs(prn, "\nStep 1 results\n");
365 pputs(prn, "unit delta s2e s2y\n");
366 }
367
368 bigrow = 0;
369
370 for (i=0; i<N && !err; i++) {
371 double yti, yti_1, delta;
372 int p_i, T_i, k_i;
373 int pt0;
374
375 if (p_varies) {
376 p_i = plist[i+1];
377 T_i = t2 - t1 + 1 - (1 + p_i);
378 k_i = m + p_i;
379 gretl_matrix_reuse(y, T_i, 1);
380 gretl_matrix_reuse(X, T_i, k_i);
381 gretl_matrix_reuse(b, k_i, 1);
382 gretl_matrix_reuse(ei, T_i, 1);
383 gretl_matrix_reuse(vi, T_i, 1);
384 } else {
385 p_i = p;
386 T_i = T;
387 k_i = k;
388 }
389
390 /* indices into Z array */
391 pt1 = t1 + (i + u0) * dset->pd;
392 pt2 = t2 + (i + u0) * dset->pd;
393 pt0 = pt1 + 1 + p_i;
394
395 /* build (full length) \delta y_t in dy */
396 s = 0;
397 for (t=pt1+1; t<=pt2; t++) {
398 yti = dset->Z[vnum][t];
399 yti_1 = dset->Z[vnum][t-1];
400 gretl_vector_set(dy, s++, yti - yti_1);
401 }
402
403 /* build y_{t-1} in y */
404 s = 0;
405 for (t=pt0; t<=pt2; t++) {
406 yti_1 = dset->Z[vnum][t-1];
407 gretl_vector_set(y, s++, yti_1);
408 }
409
410 /* augmented case: write lags of dy into X */
411 for (j=1; j<=p_i; j++) {
412 int col = m + j - 2;
413 double dylag;
414
415 s = 0;
416 for (t=pt0; t<=pt2; t++) {
417 dylag = gretl_vector_get(dy, t - pt1 - 1 - j);
418 gretl_matrix_set(X, s++, col, dylag);
419 }
420 }
421
422 /* set lagged y as last regressor */
423 for (t=0; t<T_i; t++) {
424 gretl_matrix_set(X, t, k_i - 1, y->val[t]);
425 }
426
427 #if LLC_DEBUG > 1
428 gretl_matrix_print(dy, "dy");
429 gretl_matrix_print(y, "y1");
430 gretl_matrix_print(X, "X");
431 #endif
432
433 if (p_i > 0) {
434 /* "virtual trimming" of dy for regressions
435 note: we reset these values below
436 */
437 dy->val += p_i;
438 dy->rows -= p_i;
439 }
440
441 /* run (A)DF regression: LLC eqn (1') */
442 err = gretl_matrix_ols(dy, X, b, NULL, ui, NULL);
443 if (err) {
444 break;
445 }
446
447 /* for verbose reporting */
448 delta = b->val[b->rows-1];
449
450 #if LLC_DEBUG > 1
451 gretl_matrix_print(b, "ADF coeffs");
452 #endif
453
454 if (k_i > 1) {
455 /* reduced regressor matrix for auxiliary regressions:
456 omit the last column containing the lagged level of y
457 */
458 gretl_matrix_reuse(X, T_i, k_i - 1);
459 gretl_matrix_reuse(b, k_i - 1, 1);
460
461 err = gretl_matrix_ols(dy, X, b, NULL, ei, NULL);
462 if (!err) {
463 err = gretl_matrix_ols(y, X, b, NULL, vi, NULL);
464 }
465
466 gretl_matrix_reuse(X, T, k);
467 gretl_matrix_reuse(b, k, 1);
468 } else {
469 /* no auxiliary regressions required */
470 gretl_matrix_copy_values(ei, dy);
471 gretl_matrix_copy_values(vi, y);
472 }
473
474 if (p_i > 0) {
475 /* restore dy to full length */
476 dy->val -= p_i;
477 dy->rows += p_i;
478 }
479
480 if (!err) {
481 double sui, s2yi, s2ui = 0.0;
482
483 for (t=0; t<T_i; t++) {
484 s2ui += ui->val[t] * ui->val[t];
485 }
486
487 s2ui /= T_i;
488 sui = sqrt(s2ui);
489
490 /* write normalized per-unit ei and vi into big matrices */
491 gretl_matrix_divide_by_scalar(ei, sui);
492 gretl_matrix_divide_by_scalar(vi, sui);
493 gretl_matrix_inscribe_matrix(e, ei, bigrow, 0, GRETL_MOD_NONE);
494 gretl_matrix_inscribe_matrix(v, vi, bigrow, 0, GRETL_MOD_NONE);
495 bigrow += T_i;
496
497 s2yi = LLC_lrvar(dy, K, m, &err);
498 if (!err) {
499 /* cumulate ratio of LR std dev to innovation std dev */
500 SN += sqrt(s2yi) / sui;
501 }
502
503 if (verbose) {
504 pprintf(prn, "%3d %#10.5g %#10.5g %#10.5g\n",
505 i+1, delta, s2ui, s2yi);
506 }
507 }
508
509 if (p_varies) {
510 gretl_matrix_reuse(y, T, 1);
511 gretl_matrix_reuse(X, T, k);
512 gretl_matrix_reuse(b, k, 1);
513 gretl_matrix_reuse(ei, T, 1);
514 gretl_matrix_reuse(vi, T, 1);
515 }
516 }
517
518 if (!err) {
519 /* the final step: pooled regression of e on v */
520 double delta, s2e, STD, td;
521 double mstar, sstar;
522
523 gretl_matrix_reuse(b, 1, 1);
524 err = gretl_matrix_ols(e, v, b, NULL, eps, NULL);
525
526 if (!err) {
527 double ee = 0, vv = 0;
528
529 for (t=0; t<NT; t++) {
530 ee += eps->val[t] * eps->val[t];
531 vv += v->val[t] * v->val[t];
532 }
533
534 SN /= N;
535 delta = b->val[0];
536 s2e = ee / NT;
537 STD = sqrt(s2e / vv);
538 td = delta / STD;
539
540 /* fetch the Levin-Lin-Chu correction factors */
541 err = get_LLC_corrections(T, m, &mstar, &sstar);
542 }
543
544 if (!err) {
545 double z = (td - NT * (SN / s2e) * STD * mstar) / sstar;
546 double pval = normal_cdf(z);
547
548 pprintf(prn, "\nS_N = %g, mu* = %g, s* = %g\n", SN, mstar, sstar);
549
550 if (!(opt & OPT_Q)) {
551 const char *heads[] = {
552 N_("coefficient"),
553 N_("t-ratio"),
554 N_("z-score")
555 };
556 const char *s = dset->varname[vnum];
557 char NTstr[32];
558 int sp[3] = {0, 3, 5};
559 int w[3] = {4, 6, 0};
560
561 pputc(prn, '\n');
562 pprintf(prn, _("Levin-Lin-Chu pooled ADF test for %s\n"), s);
563 pprintf(prn, "%s ", _(DF_test_spec(m)));
564
565 if (p_varies) {
566 pprintf(prn, _("including %.2f lags of (1-L)%s (average)"), pbar, s);
567 } else if (p == 1) {
568 pprintf(prn, _("including one lag of (1-L)%s"), s);
569 } else {
570 pprintf(prn, _("including %d lags of (1-L)%s"), p, s);
571 }
572 pputc(prn, '\n');
573
574 pprintf(prn, _("Bartlett truncation at %d lags\n"), K);
575 sprintf(NTstr, "N,T = (%d,%d)", N, dyT + 1);
576 pprintf(prn, _("%s, using %d observations"), NTstr, NT);
577
578 pputs(prn, "\n\n");
579 for (i=0; i<3; i++) {
580 pputs(prn, _(heads[i]));
581 bufspace(w[i], prn);
582 w[i] = sp[i] + g_utf8_strlen(_(heads[i]), -1);
583 }
584 pputc(prn, '\n');
585
586 pprintf(prn, "%*.5g %*.3f %*.6g [%.4f]\n\n",
587 w[0], delta, w[1], td, w[2], z, pval);
588 }
589
590 record_test_result(z, pval);
591 }
592 }
593
594 gretl_matrix_block_destroy(B);
595
596 return err;
597 }
598