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 /* Note, 2016-09-26: to avoid clutter and unclarity, this version of
21 gretl_midas.c defaults to "conditional OLS" for estimation of MIDAS
22 models that comprise one or more beta or normalized exponential
23 Almon (nealmon) terms. By this we mean a combination of L-BFGS-B
24 (with constraints on the beta and/or nealmon hyper-parameters) and
25 OLS: all coefficients other than the hyper-parameters are estimated
26 via OLS conditional on the "theta" vector of hyper-parameters as
27 optimized by L-BFGS-B.
28
29 In case this proves to be a bad idea, we could go back to commit
30 8aa2c9, the step before we started working towards conditional OLS.
31 */
32
33 #include "libgretl.h"
34 #include "usermat.h"
35 #include "uservar.h"
36 #include "matrix_extra.h"
37 #include "libset.h"
38 #include "gretl_bfgs.h"
39 #include "nlspec.h"
40 #include "gretl_midas.h"
41
42 #define MIDAS_DEBUG 0
43 #define FC_DEBUG 0
44
45 enum mflags {
46 M_PRELAG = 1 << 0,
47 M_AUTO = 1 << 1,
48 M_TMPLIST = 1 << 2
49 };
50
51 enum midas_methods {
52 MDS_NLS,
53 MDS_OLS,
54 MDS_BFGS,
55 MDS_GSS
56 };
57
58 /* information on an individual MIDAS term */
59
60 struct midas_term_ {
61 char lnam0[VNAMELEN]; /* name of MIDAS list on input */
62 char lname[VNAMELEN]; /* name of MIDAS list */
63 char mname[VNAMELEN]; /* name of initial theta vector */
64 gretl_matrix *theta; /* value of initial theta vector */
65 char flags; /* values from @mflags */
66 int minlag; /* minimum lag */
67 int maxlag; /* maximum lag */
68 int type; /* type of parameterization */
69 int nparm; /* number of parameters */
70 int nlags; /* number of lag terms */
71 int *laglist; /* list of lag series */
72 };
73
74 typedef struct midas_term_ midas_term;
75
76 /* information on a MIDAS model as a whole */
77
78 struct midas_info_ {
79 const int *list; /* incoming regression list */
80 int yno; /* ID number of dependent variable */
81 int nx; /* number of low-frequency regressors */
82 midas_term *mterms; /* array of MIDAS terms */
83 int nmidas; /* number of elements in @mterms */
84 int *xlist; /* list of low-frequency regressors */
85 int *seplist; /* list of coeff separator positions */
86 gchar *pnames; /* string holding parameter names */
87 int hfslopes; /* # of high-frequency slope coeffs */
88 int nalmonp; /* number of straight Almon poly terms */
89 int nbeta1; /* number of clamped beta terms */
90 int method; /* estimation method */
91 int ldepvar; /* position of lagged dependent variable? */
92 int nobs; /* number of observations used */
93 size_t colsize; /* matrix column size in bytes */
94 DATASET *dset; /* dataset */
95 gretl_matrix *b; /* all, or OLS, coefficients */
96 gretl_matrix *g; /* SSR gradients */
97 gretl_matrix *theta; /* MIDAS hyper-params, combined */
98 gretl_matrix *bounds; /* bounds on hyper-params */
99 gretl_matrix *y; /* dependent vector */
100 gretl_matrix *u; /* residual vector */
101 gretl_matrix *xi; /* single column of X */
102 gretl_matrix *X; /* X(\beta) */
103 double SSR; /* sum of squared residuals */
104 double orig_toler; /* prior NLS_TOLER value */
105 };
106
107 typedef struct midas_info_ midas_info;
108
109 #define prelag(m) (m->flags & M_PRELAG)
110 #define auto_theta(m) (m->flags & M_AUTO)
111
112 static midas_term *mterms_from_array (gretl_array *A,
113 int *nmidas,
114 int *err);
115
116 static double cond_ols_callback (double *theta, double *g,
117 int n, void *ptr);
118
119 /* Identify parameterizations which take an additional
120 leading coefficient: all but U-MIDAS and the plain
121 Almon polynomial specification.
122 */
123 #define takes_coeff(t) (t != MIDAS_U && t != MIDAS_ALMONP)
124
125 /* convenience macro for picking out beta specifications */
126 #define beta_type(t) (t == MIDAS_BETA0 || t == MIDAS_BETA1 || t == MIDAS_BETAN)
127
128 /* days per month or quarter: maybe make this user-
129 configurable? */
130
midas_days_per_period(int days_per_week,int pd)131 int midas_days_per_period (int days_per_week, int pd)
132 {
133 int ret = 0;
134
135 if (days_per_week == 5) {
136 ret = 22;
137 } else if (days_per_week == 6) {
138 ret = 26;
139 } else if (days_per_week == 7) {
140 ret = 30;
141 }
142
143 return (pd == 12)? ret : 3 * ret;
144 }
145
midas_days_to_freq(int pd,int m)146 static int midas_days_to_freq (int pd, int m)
147 {
148 int days[] = {22, 26, 30};
149 int freq[] = {5, 6, 7};
150 int i;
151
152 for (i=0; i<3; i++) {
153 if ((pd == 12 && m == days[i]) ||
154 (pd == 4 && m == 3 * days[i])) {
155 return freq[i];
156 }
157 }
158
159 return 0;
160 }
161
162 /* Could @m be a valid frequency ratio (= number of members
163 of a valid "MIDAS list"), given the periodicity of @dset?
164 If so, return 1; if not, return 0.
165 */
166
is_valid_midas_frequency_ratio(const DATASET * dset,int m)167 int is_valid_midas_frequency_ratio (const DATASET *dset, int m)
168 {
169 if (dset->pd == 1) {
170 /* lf = annual, hf = quarterly or monthly */
171 return (m == 4 || m == 12);
172 } else if (dset->pd == 4 && m == 3) {
173 /* lf = quarterly, hf = monthly */
174 return 1;
175 } else if (dset->pd == 4 || dset->pd == 12) {
176 /* lf = quarterly or monthly, hf = daily */
177 if (m == midas_days_per_period(5, dset->pd)) {
178 return 1;
179 } else if (m == midas_days_per_period(6, dset->pd)) {
180 return 1;
181 } else if (m == midas_days_per_period(7, dset->pd)) {
182 return 1;
183 }
184 }
185
186 return 0;
187 }
188
get_midas_frequency(const DATASET * dset,int m)189 int get_midas_frequency (const DATASET *dset, int m)
190 {
191 if (dset->pd == 1 && (m == 4 || m == 12)) {
192 /* lf = annual -> hf quarterly or monthly */
193 return m;
194 } else if (dset->pd == 4 && m == 3) {
195 /* lf = quarterly -> hf monthly */
196 return 12;
197 } else if (dset->pd == 4 || dset->pd == 12) {
198 /* hf daily 5, 6 or 7? */
199 return midas_days_to_freq(dset->pd, m);
200 }
201
202 return 0;
203 }
204
midas_pdstr(const DATASET * dset,int cfac)205 const char *midas_pdstr (const DATASET *dset, int cfac)
206 {
207 const char *pdstrs[] = {
208 N_("quarterly"),
209 N_("monthly"),
210 N_("daily")
211 };
212 int f = get_midas_frequency(dset, cfac);
213
214 if (f == 4) {
215 return pdstrs[0];
216 } else if (f == 12) {
217 return pdstrs[1];
218 } else if (f == 5 || f == 6 || f == 7) {
219 return pdstrs[2];
220 } else {
221 return "unknown frequency";
222 }
223 }
224
midas_m_from_pd(const DATASET * dset,int pd)225 int midas_m_from_pd (const DATASET *dset, int pd)
226 {
227 if (dset->pd == 1 && (pd == 4 || pd == 12)) {
228 return pd;
229 } else if (dset->pd == 4 && pd == 12) {
230 return 3;
231 } else if (dset->pd == 4 || dset->pd == 12) {
232 return midas_days_per_period(pd, dset->pd);
233 }
234
235 return 0;
236 }
237
238 /* Infer the current month from the current @qtr along
239 with the number of days per period, @ndays, and the
240 current index within the month-days array, @day.
241 */
242
quarter_to_month(int qtr,int ndays,int day)243 static int quarter_to_month (int qtr, int ndays, int day)
244 {
245 return qtr * 3 - 2 + (day - 1) / (ndays/3);
246 }
247
get_midas_pd(int pd,int m,int * err)248 static int get_midas_pd (int pd, int m, int *err)
249 {
250 int mpd = 0;
251
252 if (pd == 1) {
253 /* annual: midas series should be quarterly or monthly */
254 if (m != 4 && m != 12) {
255 *err = E_INVARG;
256 } else {
257 mpd = m;
258 }
259 } else if (pd == 4) {
260 /* quarterly: midas series should be monthly or daily */
261 if (m == 3) {
262 mpd = 12;
263 } else if (m == midas_days_per_period(5, 4)) {
264 mpd = 5;
265 } else if (m == midas_days_per_period(6, 4)) {
266 mpd = 6;
267 } else if (m == midas_days_per_period(7, 4)) {
268 mpd = 7;
269 } else {
270 *err = E_INVARG;
271 }
272 } else {
273 /* monthly: midas series should be daily */
274 if (m == midas_days_per_period(5, 12)) {
275 mpd = 5;
276 } else if (m == midas_days_per_period(6, 12)) {
277 mpd = 6;
278 } else if (m == midas_days_per_period(7, 12)) {
279 mpd = 7;
280 } else {
281 *err = E_INVARG;
282 }
283 }
284
285 return mpd;
286 }
287
288 /* Construct an auxiliary dataset in which the data from
289 a MIDAS list are represented at their "native"
290 frequency. We use this for pretty-printing a high
291 frequency series.
292 */
293
midas_aux_dataset(const int * list,const DATASET * dset,int * err)294 DATASET *midas_aux_dataset (const int *list,
295 const DATASET *dset,
296 int *err)
297 {
298 DATASET *mset = NULL;
299 gretlopt opt = 0;
300 int mpd, pd = dset->pd;
301 int T, m = list[0];
302 int yr, mon;
303 int daily = 0;
304
305 if (m < 3 || gretl_list_has_separator(list)) {
306 *err = E_INVARG;
307 } else if (!dataset_is_time_series(dset)) {
308 *err = E_INVARG;
309 } else if (pd != 1 && pd != 4 && pd != 12) {
310 /* host dataset should be annual, quarterly or monthly */
311 *err = E_PDWRONG;
312 }
313
314 if (!*err) {
315 mpd = get_midas_pd(pd, m, err);
316 }
317
318 if (*err) {
319 return NULL;
320 }
321
322 if (!gretl_is_midas_list(list, dset)) {
323 gretl_warnmsg_set(_("The argument does not seem to be a MIDAS list"));
324 }
325
326 T = sample_size(dset) * m;
327
328 if (mpd >= 5 && mpd <= 7) {
329 /* we'll add markers for daily dates */
330 daily = 1;
331 opt = OPT_M;
332 }
333
334 mset = create_auxiliary_dataset(1, T, opt);
335 if (mset == NULL) {
336 *err = E_ALLOC;
337 }
338
339 if (!*err) {
340 char *p, obs[OBSLEN];
341 int nonex, qtr = 0;
342 int i, t, s, m3 = 0;
343
344 mset->pd = mpd;
345 mset->structure = TIME_SERIES;
346 strcpy(mset->varname[0], dset->varname[list[1]]);
347 p = strrchr(mset->varname[0], '_');
348 if (p != NULL) *p = '\0';
349
350 ntolabel(obs, dset->t1, dset);
351
352 if (mpd == 4) {
353 sprintf(mset->stobs, "%d:1", atoi(obs));
354 } else if (mpd == 12) {
355 sprintf(mset->stobs, "%d:01", atoi(obs));
356 }
357
358 if (daily && pd == 4) {
359 m3 = m / 3;
360 }
361
362 /* loop across observations in low-frequency dataset */
363
364 s = 0;
365 for (t=dset->t1; t<=dset->t2; t++) {
366 if (daily) {
367 ntolabel(obs, t, dset);
368 sscanf(obs, "%d:%d", &yr, &mon);
369 if (pd == 4) {
370 qtr = mon;
371 }
372 }
373 /* read data right-to-left! */
374 for (i=m; i>0; i--) {
375 int vi = list[i];
376
377 if (daily) {
378 if (pd == 4) {
379 mon = quarter_to_month(qtr, m, m-i+1);
380 nonex = daily_index_to_date(mset->S[s], yr, mon,
381 (m-i) % m3, mpd);
382 } else {
383 nonex = daily_index_to_date(mset->S[s], yr, mon,
384 m-i, mpd);
385 }
386 if (nonex) {
387 /* skip any non-existent daily dates */
388 mset->t2 -= 1;
389 } else {
390 mset->Z[0][s++] = dset->Z[vi][t];
391 }
392 } else {
393 mset->Z[0][s++] = dset->Z[vi][t];
394 }
395 }
396 }
397
398 if (daily) {
399 strcpy(mset->stobs, mset->S[0]);
400 strcpy(mset->endobs, mset->S[mset->t2]);
401 mset->markers = DAILY_DATE_STRINGS;
402 }
403
404 mset->sd0 = get_date_x(mset->pd, mset->stobs);
405 if (!daily) {
406 ntolabel(mset->endobs, mset->t2, mset);
407 }
408 }
409
410 return mset;
411 }
412
midas_list_check(const int * list,const DATASET * dset)413 static int midas_list_check (const int *list,
414 const DATASET *dset)
415 {
416 int mpd, m = list[0];
417 int err = 0;
418
419 if (m < 3 || gretl_list_has_separator(list)) {
420 return E_INVARG;
421 } else if (!dataset_is_time_series(dset)) {
422 return E_PDWRONG;
423 } else if (dset->pd != 1 && dset->pd != 4) {
424 /* host dataset should be annual or quarterly */
425 return E_PDWRONG;
426 }
427
428 mpd = get_midas_pd(dset->pd, m, &err);
429 if (!err && mpd != 3 && mpd != 4 && mpd != 12) {
430 err = E_INVARG;
431 }
432
433 return err;
434 }
435
midas_list_to_vector(const int * list,const DATASET * dset,int * err)436 gretl_matrix *midas_list_to_vector (const int *list,
437 const DATASET *dset,
438 int *err)
439 {
440 gretl_matrix *mvec;
441 int t, i, vi, n_tail = 0;
442 int T, m = list[0];
443
444 *err = midas_list_check(list, dset);
445 if (*err) {
446 return NULL;
447 }
448
449 if (dset->t2 < dset->n - 1) {
450 t = dset->t2 + 1;
451 for (i=m; i>0; i--) {
452 vi = list[i];
453 if (!na(dset->Z[vi][t])) {
454 n_tail++;
455 } else {
456 break;
457 }
458 }
459 }
460
461 T = sample_size(dset) * m + n_tail;
462 mvec = gretl_matrix_alloc(T, 1);
463
464 if (mvec == NULL) {
465 *err = E_ALLOC;
466 } else {
467 int s = 0;
468
469 for (t=dset->t1; t<=dset->t2; t++) {
470 /* read data right-to-left! */
471 for (i=m; i>0; i--) {
472 vi = list[i];
473 mvec->val[s++] = dset->Z[vi][t];
474 }
475 }
476 if (n_tail > 0) {
477 t = dset->t2 + 1;
478 for (i=0; i<n_tail; i++) {
479 vi = list[m-i];
480 mvec->val[s++] = dset->Z[vi][t];
481 }
482 }
483 }
484
485 return mvec;
486 }
487
make_auto_theta(char * mstr,int i,int ptype,int k,int m1,int m2,int * err)488 static gretl_matrix *make_auto_theta (char *mstr, int i,
489 int ptype, int k,
490 int m1, int m2,
491 int *err)
492 {
493 gretl_matrix *theta = NULL;
494
495 if (k == 0) {
496 /* we have to infer k */
497 if (*mstr == '\0' || !strcmp(mstr, "null")) {
498 /* OK if we know how many parameters are needed? */
499 if (ptype == MIDAS_BETA0 || ptype == MIDAS_BETA1) {
500 k = 2;
501 } else if (ptype == MIDAS_BETAN) {
502 k = 3;
503 } else if (ptype == MIDAS_NEALMON) {
504 /* we'll default to 2 */
505 k = 2;
506 } else if (ptype == MIDAS_U) {
507 k = m2 - m1 + 1;
508 } else {
509 *err = E_INVARG;
510 }
511 } else if (integer_string(mstr)) {
512 /* does @mstr give the number of params? */
513 int chk = atoi(mstr);
514
515 if (chk >= 1 && chk < 100) {
516 k = chk;
517 }
518 } else {
519 /* try treating @mstr as matrix definition? */
520 theta = generate_matrix(mstr, NULL, err);
521 }
522 }
523
524 if (k > 0) {
525 theta = gretl_zero_matrix_new(1, k);
526 if (theta != NULL) {
527 if (beta_type(ptype)) {
528 theta->val[0] = 1;
529 theta->val[1] = 10; /* optimize somehow? */
530 }
531 } else {
532 *err = E_ALLOC;
533 }
534 }
535
536 if (theta != NULL) {
537 sprintf(mstr, "theta___%d", i+1);
538 private_matrix_add(theta, mstr);
539 } else if (!*err) {
540 *err = E_PARSE;
541 }
542
543 #if MIDAS_DEBUG
544 gretl_matrix_print(theta, "auto-generated theta");
545 #endif
546
547 return theta;
548 }
549
lag_info_from_prelag_list(midas_term * mt,const int * list,const DATASET * dset)550 static int lag_info_from_prelag_list (midas_term *mt,
551 const int *list,
552 const DATASET *dset)
553 {
554 int mf = series_get_midas_freq(dset, list[1]);
555 int m1 = series_get_lag(dset, list[1]);
556 int p1 = series_get_midas_period(dset, list[1]);
557 int i, p, maxp = 0;
558
559 if (mf > 0) {
560 maxp = midas_m_from_pd(dset, mf);
561 }
562
563 if (maxp == 0 && p1 > 0) {
564 for (i=1; i<=list[0]; i++) {
565 p = series_get_midas_period(dset, list[i]);
566 if (p > maxp) {
567 maxp = p;
568 }
569 }
570 }
571
572 if (is_valid_midas_frequency_ratio(dset, maxp)) {
573 int hfl = maxp * m1 - (p1 - 1);
574
575 mt->minlag = hfl;
576 mt->maxlag = hfl + list[0] - 1;
577 } else {
578 /* oof! just report 1,... */
579 mt->minlag = 1;
580 mt->maxlag = list[0];
581 }
582
583 return 0;
584 }
585
midas_term_init(midas_term * mt)586 static void midas_term_init (midas_term *mt)
587 {
588 mt->lnam0[0] = '\0';
589 mt->lname[0] = '\0';
590 mt->mname[0] = '\0';
591 mt->theta = NULL;
592 mt->flags = 0;
593 mt->minlag = 0;
594 mt->maxlag = 0;
595 mt->type = 0;
596 mt->nparm = 0;
597 mt->nlags = 0;
598 }
599
get_p1_p2(const char * s1,const char * s2,int * p1,int * p2)600 static int get_p1_p2 (const char *s1, const char *s2,
601 int *p1, int *p2)
602 {
603 int err = 0;
604
605 *p1 = gretl_int_from_string(s1, &err);
606
607 if (!err) {
608 *p2 = gretl_int_from_string(s2, &err);
609 }
610
611 return err;
612 }
613
midas_type_from_string(const char * s)614 static int midas_type_from_string (const char *s)
615 {
616 if (!strcmp(s, "\"umidas\"") || !strcmp(s, "umidas")) {
617 return MIDAS_U;
618 } else if (!strcmp(s, "\"nealmon\"") || !strcmp(s, "nealmon")) {
619 return MIDAS_NEALMON;
620 } else if (!strcmp(s, "\"beta0\"") || !strcmp(s, "beta0")) {
621 return MIDAS_BETA0;
622 } else if (!strcmp(s, "\"betan\"") || !strcmp(s, "betan")) {
623 return MIDAS_BETAN;
624 } else if (!strcmp(s, "\"almonp\"") || !strcmp(s, "almonp")) {
625 return MIDAS_ALMONP;
626 } else if (!strcmp(s, "\"beta1\"") || !strcmp(s, "beta1")) {
627 return MIDAS_BETA1;
628 } else {
629 return -1;
630 }
631 }
632
midas_type_from_arg(const char * s,int * umidas,int * err)633 static int midas_type_from_arg (const char *s, int *umidas, int *err)
634 {
635 int ret = -1;
636
637 if (integer_string(s)) {
638 ret = atoi(s);
639 } else {
640 ret = midas_type_from_string(s);
641 if (ret < 0) {
642 char *sval = get_string_by_name(s);
643
644 if (sval != NULL) {
645 ret = midas_type_from_string(sval);
646 }
647 }
648 if (ret < 0 && gretl_is_scalar(s)) {
649 ret = gretl_scalar_get_value(s, err);
650 }
651 }
652
653 if (ret < MIDAS_U || ret >= MIDAS_MAX) {
654 *err = E_INVARG;
655 } else if (ret == 0) {
656 *umidas = 1;
657 }
658
659 return ret;
660 }
661
662 /* Parse a particular entry in the incoming array of MIDAS
663 specifications. Each entry should look like one of
664 the following:
665
666 (1) mds(list, minlag, maxlag, type, theta)
667 (2) mds(list, minlag, maxlag, 0)
668
669 (3) mdsl(list, type, theta)
670 (4) mdsl(list, 0)
671
672 Forms (1) and (2) imply that lags of the MIDAS terms
673 in @list should be auto-generated, while (3) and (4)
674 imply that the incoming @list already holds whatever
675 lags are wanted.
676
677 Variants (2) and (4), with @type set to 0, are for
678 U-MIDAS terms: in this case @theta (a vector to
679 initialize the coefficients) is not required.
680 */
681
parse_midas_term(const char * s,midas_term * mt,int i,const DATASET * dset)682 static int parse_midas_term (const char *s,
683 midas_term *mt,
684 int i,
685 const DATASET *dset)
686 {
687 char lname[VNAMELEN];
688 char mname[VNAMELEN];
689 char p1str[VNAMELEN];
690 char p2str[VNAMELEN];
691 char typestr[10];
692 char fmt[64];
693 int ns, p1 = 0, p2 = 0;
694 int type, umidas = 0;
695 int err = 0;
696
697 midas_term_init(mt);
698 *typestr = *lname = *mname = '\0';
699
700 if (!strncmp(s, "mds(", 4)) {
701 /* calling for auto-generated lags */
702 s += 4;
703 sprintf(fmt, "%%%d[^, ] , %%%d[^, ] , %%%d[^, ] , %%%d[^,) ] , %%%d[^) ])",
704 VNAMELEN-1, VNAMELEN-1, VNAMELEN-1, 9, VNAMELEN-1);
705 ns = sscanf(s, fmt, lname, p1str, p2str, typestr, mname);
706 err = get_p1_p2(p1str, p2str, &p1, &p2);
707 if (!err) {
708 type = midas_type_from_arg(typestr, &umidas, &err);
709 }
710 if (!err && ns < 4) {
711 err = E_PARSE;
712 }
713 } else if (!strncmp(s, "mdsl(", 5)) {
714 /* list already holds lags */
715 mt->flags |= M_PRELAG;
716 s += 5;
717 sprintf(fmt, "%%%d[^, ] , %%%d[^,) ] , %%%d[^) ])",
718 VNAMELEN-1, 9, VNAMELEN-1);
719 ns = sscanf(s, fmt, lname, typestr, mname);
720 type = midas_type_from_arg(typestr, &umidas, &err);
721 if (!err && ns < 2) {
722 err = E_PARSE;
723 }
724 } else {
725 err = E_INVARG;
726 }
727
728 if (!err) {
729 int *list = get_list_by_name(lname);
730 int k = 0;
731
732 if (!umidas) {
733 gretl_matrix *theta = NULL;
734 char c = mname[0];
735
736 if (c != '\0' && !isdigit(c) && c != '{' && strcmp(mname, "null")) {
737 theta = get_matrix_by_name(mname);
738 }
739
740 if (theta != NULL) {
741 /* copy, don't clobber, the user's initializer */
742 sprintf(mname, "theta___%d", i+1);
743 mt->theta = gretl_matrix_copy(theta);
744 if (mt->theta != NULL) {
745 private_matrix_add(mt->theta, mname);
746 }
747 } else {
748 /* create automatic initializer */
749 mt->flags |= M_AUTO;
750 mt->theta = make_auto_theta(mname, i, type, 0,
751 p1, p2, &err);
752 }
753 }
754
755 if (err) {
756 ; /* leave it alone */
757 } else if (!umidas && mt->theta == NULL) {
758 err = E_ALLOC;
759 } else if (prelag(mt) && list == NULL) {
760 err = E_INVARG;
761 } else if (!prelag(mt) && !gretl_is_midas_list(list, dset)) {
762 gretl_errmsg_set("mds(): the first term must be a MIDAS list");
763 err = E_INVARG;
764 } else if (p1 > p2) {
765 err = E_INVARG;
766 } else if (type < 0 || type >= MIDAS_MAX) {
767 err = E_INVARG;
768 } else if (umidas) {
769 if (prelag(mt)) {
770 k = list[0];
771 } else {
772 k = p2 - p1 + 1;
773 }
774 } else {
775 k = gretl_vector_get_length(mt->theta);
776 if (k < 1 || (type == MIDAS_BETA0 && k != 2) ||
777 (type == MIDAS_BETA1 && k != 2) ||
778 (type == MIDAS_BETAN && k != 3)) {
779 err = E_INVARG;
780 }
781 }
782
783 if (!err) {
784 strcpy(mt->lnam0, lname);
785 strcpy(mt->lname, lname);
786 if (!umidas) {
787 strcpy(mt->mname, mname);
788 }
789 if (prelag(mt)) {
790 /* scrounge lag info from incoming list */
791 lag_info_from_prelag_list(mt, list, dset);
792 } else {
793 mt->minlag = p1;
794 mt->maxlag = p2;
795 }
796 mt->type = type;
797 mt->nparm = k;
798 }
799 }
800
801 return err;
802 }
803
804 /* In case we got any U-MIDAS terms in the specification,
805 check to see if we need to add any initializers for
806 the associated coefficients.
807 */
808
umidas_check(midas_info * mi,int n_umidas)809 static int umidas_check (midas_info *mi, int n_umidas)
810 {
811 int i, err = 0;
812
813 if (n_umidas == mi->nmidas) {
814 /* all U-MIDAS: so use OLS */
815 mi->method = MDS_OLS;
816 return 0;
817 }
818
819 /* mix of U-MIDAS and "unboxed" terms: we'll be using
820 NLS and we need initializers */
821
822 for (i=0; i<mi->nmidas && !err; i++) {
823 midas_term *mt = &mi->mterms[i];
824
825 if (mt->type == MIDAS_U && mt->theta == NULL) {
826 mt->flags |= M_AUTO;
827 mt->theta = make_auto_theta(mt->mname, i, MIDAS_U,
828 mt->nparm, 0, 0, &err);
829 }
830 }
831
832 return err;
833 }
834
835 /* Parse the @spec string, which should contain one or more
836 MIDAS specifications. For details on what exactly we're
837 looking for, see the comment on parse_midas_term() above.
838 This function also potentially revises the estimation
839 method: to OLS if we have nothing but UMIDAS terms; to
840 MDS_BFGS if we "boxed" terms and the --levenberg option
841 has not been specified.
842 */
843
844 static int
parse_midas_specs(midas_info * mi,const char * spec,const DATASET * dset,gretlopt opt)845 parse_midas_specs (midas_info *mi, const char *spec,
846 const DATASET *dset, gretlopt opt)
847 {
848 const char *s;
849 int n_spec = 0;
850 int n_umidas = 0;
851 int n_boxed = 0;
852 int n_beta1 = 0;
853 int n_almonp = 0;
854 int err = 0;
855
856 /* first check: count closing parentheses */
857 s = spec;
858 while (*s) {
859 if (*s == ')') {
860 n_spec++;
861 }
862 s++;
863 }
864
865 if (n_spec == 0) {
866 /* spec is junk! */
867 err = E_PARSE;
868 } else {
869 /* allocate info structs */
870 mi->mterms = malloc(n_spec * sizeof *mi->mterms);
871 if (mi->mterms == NULL) {
872 err = E_ALLOC;
873 }
874 }
875
876 /* number of slope coeffs on high-frequency terms */
877 mi->hfslopes = 0;
878
879 if (!err) {
880 /* parse and record individual MIDAS specs */
881 char test[128];
882 const char *p;
883 int len, i = 0;
884
885 s = spec;
886 for (i=0; i<n_spec && !err; i++) {
887 midas_term *mt = &mi->mterms[i];
888
889 while (*s == ' ') s++;
890 p = s;
891 while (*p && *p != ')') p++;
892 len = p - s + 1;
893 if (len > 127) {
894 err = E_PARSE;
895 } else {
896 *test = '\0';
897 strncat(test, s, len);
898 err = parse_midas_term(test, mt, i, dset);
899 if (!err) {
900 if (mt->type == MIDAS_BETA0 && (opt & OPT_C)) {
901 /* support legacy option */
902 mt->type = MIDAS_BETA1;
903 }
904 if (mt->type == MIDAS_U) {
905 n_umidas++;
906 } else if (beta_type(mt->type) ||
907 mt->type == MIDAS_NEALMON) {
908 n_boxed++;
909 } else if (mt->type == MIDAS_ALMONP) {
910 n_almonp++;
911 }
912 if (takes_coeff(mt->type)) {
913 mi->hfslopes += 1;
914 }
915 if (mt->type == MIDAS_BETA1) {
916 n_beta1++;
917 }
918 }
919 s = p + 1;
920 }
921 }
922 }
923
924 if (!err && n_beta1 > 0) {
925 /* check validity of beta clamping */
926 if (n_almonp > 0) {
927 gretl_errmsg_set("Can't combine beta1 with almonp terms");
928 err = E_INVARG;
929 } else if (opt & OPT_L) {
930 /* can't combine beta1 and --levenberg */
931 err = E_BADOPT;
932 }
933 }
934
935 if (!err) {
936 mi->nmidas = n_spec;
937 mi->nalmonp = n_almonp;
938 mi->nbeta1 = n_beta1;
939 if (n_boxed > 0 && n_almonp == 0) {
940 /* prefer LBFGS, but respect OPT_L if it's given */
941 if (!(opt & OPT_L)) {
942 /* ! --levenberg */
943 mi->method = MDS_BFGS;
944 }
945 } else if (n_umidas > 0) {
946 /* note: switch to MDS_OLS if appropriate */
947 err = umidas_check(mi, n_umidas);
948 }
949 }
950
951 if (err) {
952 /* scrub all this stuff to avoid a crash on cleanup */
953 free(mi->mterms);
954 mi->mterms = NULL;
955 mi->nmidas = 0;
956 mi->nalmonp = 0;
957 mi->nbeta1 = 0;
958 }
959
960 return err;
961 }
962
963 /* Extract the list of regular (low-frequency) regressors.
964 While we're at it, see if any regressors are lags of the
965 dependent variable and if so record this fact in the
966 ldepvar member of @mi.
967 */
968
make_midas_xlist(midas_info * mi)969 static int make_midas_xlist (midas_info *mi)
970 {
971 int ldvpos = 0, ldv_count = 0;
972 int i, xi, err = 0;
973
974 mi->nx = mi->list[0] - 1;
975
976 if (mi->nx > 0) {
977 mi->xlist = gretl_list_new(mi->nx);
978 if (mi->xlist == NULL) {
979 err = E_ALLOC;
980 } else {
981 for (i=2; i<=mi->list[0]; i++) {
982 xi = mi->list[i];
983 mi->xlist[i-1] = xi;
984 if (standard_lag_of(xi, mi->yno, mi->dset) > 0) {
985 ldvpos = i;
986 ldv_count++;
987 }
988 }
989 }
990 }
991
992 if (ldv_count == 1) {
993 mi->ldepvar = ldvpos;
994 }
995
996 return err;
997 }
998
999 /* Given the name of an incoming MIDAS list plus minlag and
1000 maxlag values (at this point stored in the midas_term
1001 structure, @mt) build the list of required lags of the
1002 MIDAS series. Or in case the incoming list already
1003 includes the required lags, just take a pointer to it.
1004 */
1005
make_midas_laglist(midas_term * mt,DATASET * dset,int * err)1006 static int *make_midas_laglist (midas_term *mt,
1007 DATASET *dset,
1008 int *err)
1009 {
1010 int *list = get_list_by_name(mt->lname);
1011
1012 if (list == NULL) {
1013 fprintf(stderr, "make_midas_laglist, '%s': no list!\n",
1014 mt->lname);
1015 *err = E_DATA;
1016 return NULL;
1017 }
1018
1019 if (!prelag(mt) && mt->minlag == 1 && mt->maxlag == list[0]) {
1020 /* the incoming list was not flagged as "mdsl", and so
1021 must be a MIDAS list, but it nonetheless contains all
1022 the lags we need
1023 */
1024 mt->flags |= M_PRELAG;
1025 }
1026
1027 if (prelag(mt)) {
1028 /* don't copy the list (and don't free it either!) */
1029 return list;
1030 } else {
1031 /* copy, because we're going to modify the list */
1032 int *lcpy = gretl_list_copy(list);
1033
1034 if (lcpy == NULL) {
1035 *err = E_ALLOC;
1036 } else {
1037 *err = list_laggenr(&lcpy, mt->minlag, mt->maxlag,
1038 NULL, dset, lcpy[0], OPT_L);
1039 }
1040
1041 return lcpy;
1042 }
1043 }
1044
fcast_make_midas_laglists(midas_term * mterms,int nmidas,DATASET * dset)1045 static int fcast_make_midas_laglists (midas_term *mterms,
1046 int nmidas,
1047 DATASET *dset)
1048 {
1049 midas_term *mt;
1050 int *mlist;
1051 int i, err = 0;
1052
1053 for (i=0; i<nmidas && !err; i++) {
1054 mt = &mterms[i];
1055 mlist = make_midas_laglist(mt, dset, &err);
1056 if (!err) {
1057 if (!prelag(mt)) {
1058 sprintf(mt->lname, "ML___%d", i+1);
1059 err = remember_list(mlist, mt->lname, NULL);
1060 mt->flags |= M_TMPLIST;
1061 }
1062 mt->nlags = mlist[0];
1063 mt->laglist = mlist;
1064 }
1065 }
1066
1067 return err;
1068 }
1069
make_midas_laglists(midas_info * mi,DATASET * dset)1070 static int make_midas_laglists (midas_info *mi,
1071 DATASET *dset)
1072 {
1073 midas_term *mt;
1074 int *mlist;
1075 int i, err = 0;
1076
1077 for (i=0; i<mi->nmidas && !err; i++) {
1078 mt = &mi->mterms[i];
1079 mlist = make_midas_laglist(mt, dset, &err);
1080 if (!err) {
1081 /* In the "prelag" case the laglist is already
1082 in userspace. If method is BFGS the list
1083 doesn't have to be pushed into userspace,
1084 but we'll record it nonetheless for model
1085 printout purposes.
1086 */
1087 if (!prelag(mt)) {
1088 sprintf(mt->lname, "ML___%d", i+1);
1089 /* note: remember_list copies its first arg */
1090 err = remember_list(mlist, mt->lname, NULL);
1091 mt->flags |= M_TMPLIST;
1092 }
1093 mt->nlags = mlist[0];
1094 mt->laglist = mlist;
1095 }
1096 }
1097
1098 return err;
1099 }
1100
1101 /* If @list is non-NULL, build a full list of all series
1102 involved: dependent variable, regular regressors, and
1103 all lags of MIDAS terms. We want this either for setting
1104 the usable sample range, or (in the case of U-MIDAS via
1105 OLS), as the list to pass to lsq().
1106
1107 If @list is NULL, however, the returned list will just
1108 contain all the MIDAS lag terms: we use this variant
1109 when setting up for forecasting.
1110 */
1111
make_midas_biglist(const int * list,midas_term * mterms,int nmidas)1112 static int *make_midas_biglist (const int *list,
1113 midas_term *mterms,
1114 int nmidas)
1115 {
1116 int i, j, nt = 0;
1117 int *biglist = NULL;
1118
1119 if (list != NULL) {
1120 nt = list[0];
1121 }
1122
1123 for (j=0; j<nmidas; j++) {
1124 nt += mterms[j].nlags;
1125 }
1126
1127 biglist = gretl_list_new(nt);
1128
1129 if (biglist != NULL) {
1130 int k = 1;
1131
1132 if (list != NULL) {
1133 for (i=1; i<=list[0]; i++) {
1134 biglist[k++] = list[i];
1135 }
1136 }
1137 for (j=0; j<nmidas; j++) {
1138 for (i=1; i<=mterms[j].nlags; i++) {
1139 biglist[k++] = mterms[j].laglist[i];
1140 }
1141 }
1142 }
1143
1144 return biglist;
1145 }
1146
1147 /* If we're doing MIDAS via nls, set the usable
1148 sample range first. That way we'll know how
1149 big the X data matrix ought to be.
1150 */
1151
midas_set_sample(midas_info * mi,DATASET * dset)1152 static int midas_set_sample (midas_info *mi,
1153 DATASET *dset)
1154 {
1155 int *biglist;
1156 int err = 0;
1157
1158 biglist = make_midas_biglist(mi->list, mi->mterms, mi->nmidas);
1159 if (biglist == NULL) {
1160 err = E_ALLOC;
1161 }
1162
1163 if (!err) {
1164 int t1 = dset->t1, t2 = dset->t2;
1165
1166 err = list_adjust_sample(biglist, &t1, &t2, dset, NULL);
1167 if (!err) {
1168 dset->t1 = t1;
1169 dset->t2 = t2;
1170 mi->nobs = t2 - t1 + 1;
1171 mi->colsize = mi->nobs * sizeof(double);
1172 }
1173 }
1174
1175 free(biglist);
1176
1177 #if MIDAS_DEBUG
1178 fprintf(stderr, "midas_set_sample: returning %d\n", err);
1179 #endif
1180
1181 return err;
1182 }
1183
1184 /* estimate unrestricted MIDAS via OLS */
1185
umidas_ols(MODEL * pmod,midas_info * mi,DATASET * dset,gretlopt opt)1186 static int umidas_ols (MODEL *pmod,
1187 midas_info *mi,
1188 DATASET *dset,
1189 gretlopt opt)
1190 {
1191 int *biglist;
1192
1193 biglist = make_midas_biglist(mi->list, mi->mterms, mi->nmidas);
1194
1195 if (biglist == NULL) {
1196 return E_ALLOC;
1197 } else {
1198 *pmod = lsq(biglist, dset, OLS, opt | OPT_Z);
1199 free(biglist);
1200 }
1201
1202 return pmod->errcode;
1203 }
1204
1205 /* For forecasting, put into uservar space a "private"
1206 matrix containing all the coefficients on MIDAS
1207 lag terms.
1208 */
1209
push_midas_coeff_array(const MODEL * pmod)1210 static int push_midas_coeff_array (const MODEL *pmod)
1211 {
1212 gretl_matrix *mhfb, *hfb;
1213 int err = 0;
1214
1215 mhfb = gretl_model_get_data(pmod, "hfb");
1216 if (mhfb == NULL) {
1217 fprintf(stderr, "Couldn't retrieve gross midas coeffs\n");
1218 err = E_DATA;
1219 } else {
1220 hfb = gretl_matrix_copy(mhfb);
1221 if (hfb == NULL) {
1222 err = E_ALLOC;
1223 }
1224 }
1225
1226 #if FC_DEBUG
1227 gretl_matrix_print(hfb, "hfb in fcast");
1228 #endif
1229
1230 if (!err) {
1231 err = private_matrix_add(hfb, "hfb___");
1232 }
1233
1234 return err;
1235 }
1236
1237 /* Prepare for generating a MIDAS forecast: this
1238 is called from midas_fcast() in forecast.c.
1239 */
1240
midas_forecast_setup(const MODEL * pmod,DATASET * dset,ForecastMethod method,char ** pformula)1241 int midas_forecast_setup (const MODEL *pmod,
1242 DATASET *dset,
1243 ForecastMethod method,
1244 char **pformula)
1245 {
1246 gretl_array *mA;
1247 midas_term *mterms = NULL;
1248 int *xlist = NULL;
1249 int *hflist = NULL;
1250 int nmidas = 0;
1251 int nx = 0;
1252 int err = 0;
1253
1254 mA = gretl_model_get_data(pmod, "midas_info");
1255 if (mA == NULL) {
1256 err = E_DATA;
1257 }
1258
1259 if (!err && gretl_model_get_int(pmod, "no_lfx") == 0) {
1260 xlist = gretl_model_get_list(pmod, "lfxlist");
1261 if (xlist == NULL) {
1262 err = E_DATA;
1263 } else {
1264 nx = xlist[0];
1265 }
1266 }
1267
1268 if (!err) {
1269 mterms = mterms_from_array(mA, &nmidas, &err);
1270 }
1271
1272 if (!err) {
1273 /* reconstitute MIDAS lag-lists */
1274 err = fcast_make_midas_laglists(mterms, nmidas, dset);
1275 }
1276
1277 if (!err) {
1278 /* build and push list of all MIDAS terms */
1279 hflist = make_midas_biglist(NULL, mterms, nmidas);
1280 if (hflist == NULL) {
1281 err = E_ALLOC;
1282 } else {
1283 /* note: remember_list() copies its argument */
1284 #if FC_DEBUG
1285 int j;
1286 printlist(hflist, "hflist");
1287 for (j=1; j<=hflist[0]; j++) {
1288 fprintf(stderr, "%d %s\n", hflist[j], dset->varname[hflist[j]]);
1289 }
1290 #endif
1291 err = remember_list(hflist, "HFL___", NULL);
1292 user_var_privatize_by_name("HFL___");
1293 free(hflist);
1294 }
1295 }
1296
1297 if (!err) {
1298 /* retrieve and push vector of all MIDAS coeffs */
1299 err = push_midas_coeff_array(pmod);
1300 }
1301
1302 if (!err) {
1303 /* build a string that can be passed to "genr" to
1304 calculate fitted values */
1305 char tmp[64], line[MAXLEN];
1306 double *b = pmod->coeff;
1307 int p, yno = pmod->list[1];
1308 int i, xi, j = 0;
1309
1310 *line = '\0';
1311 gretl_push_c_numeric_locale();
1312
1313 for (i=0; i<nx; i++) {
1314 xi = xlist[i+1];
1315 if (xi == 0) {
1316 sprintf(tmp, "%+.15g", b[j++]);
1317 } else {
1318 /* allow for dynamic formula? */
1319 p = (method == FC_STATIC)? 0 :
1320 standard_lag_of(xi, yno, dset);
1321 if (p > 0) {
1322 sprintf(tmp, "%+.15g*%s(-%d)", b[j++],
1323 dset->varname[yno], p);
1324 } else {
1325 sprintf(tmp, "%+.15g*%s", b[j++],
1326 dset->varname[xi]);
1327 }
1328 }
1329 strcat(line, tmp);
1330 }
1331 strcat(line, "+lincomb(HFL___,hfb___)");
1332
1333 gretl_pop_c_numeric_locale();
1334
1335 #if FC_DEBUG
1336 fprintf(stderr, "formula='%s'\n", line);
1337 #endif
1338 *pformula = gretl_strdup(line);
1339 }
1340
1341 free(mterms);
1342
1343 return err;
1344 }
1345
midas_beta_init(midas_info * mi)1346 static int midas_beta_init (midas_info *mi)
1347 {
1348 midas_term *mt = &mi->mterms[0];
1349 int i;
1350
1351 if (!auto_theta(mt)) {
1352 /* the user provided a theta value, so use it */
1353 for (i=0; i<mt->nparm; i++) {
1354 mi->theta->val[i] = mt->theta->val[i];
1355 }
1356 } else {
1357 /* Trial some alternative initializations for the
1358 beta hyper-parameters (theta) and select the one
1359 that produces the smallest SSR on running OLS.
1360 */
1361 double theta[3][2] = {
1362 {1, 1}, {1, 10}, {2, 10}
1363 };
1364 double SSR, SSRmin = 1e200;
1365 int imax = mt->type == MIDAS_BETA1 ? 2 : 3;
1366 int best_idx = -1;
1367
1368 for (i=0; i<imax; i++) {
1369 SSR = cond_ols_callback(theta[i], NULL, 2, mi);
1370 if (SSR < SSRmin) {
1371 SSRmin = SSR;
1372 best_idx = i;
1373 }
1374 }
1375
1376 if (best_idx >= 0) {
1377 mi->theta->val[0] = theta[best_idx][0];
1378 mi->theta->val[1] = theta[best_idx][1];
1379 #if MIDAS_DEBUG
1380 fprintf(stderr, "midas_beta_init: best_idx = %d, SSRmin=%g\n",
1381 best_idx, SSRmin);
1382 gretl_matrix_print(mi->theta, "best theta");
1383 #endif
1384 }
1385 }
1386
1387 return 0;
1388 }
1389
1390 /* L-BFGS-B apparatus */
1391
midas_bfgs_setup(midas_info * mi,DATASET * dset,gretlopt opt)1392 static int midas_bfgs_setup (midas_info *mi, DATASET *dset,
1393 gretlopt opt)
1394 {
1395 double eps = pow(2.0, -52);
1396 double *src, *targ = NULL;
1397 int i, j, k, ii, vi;
1398 int nb, bound_rows = 0;
1399 int type0 = 0;
1400 midas_term *mt;
1401
1402 if (mi->nbeta1 == 1 && mi->nmidas == 1) {
1403 mi->method = MDS_GSS;
1404 }
1405
1406 mi->u = gretl_column_vector_alloc(mi->nobs);
1407 mi->y = gretl_column_vector_alloc(mi->nobs);
1408 mi->xi = gretl_column_vector_alloc(mi->nobs);
1409
1410 if (mi->u == NULL || mi->y == NULL || mi->xi == NULL) {
1411 return E_ALLOC;
1412 }
1413
1414 if (mi->nalmonp == 0) {
1415 memcpy(mi->y->val, dset->Z[mi->yno] + dset->t1, mi->colsize);
1416 }
1417
1418 /* number of low-frequency regressors */
1419 nb = mi->nx;
1420
1421 /* How many boxed coeffs do we have? And how many
1422 coeffs in conditional OLS (@nb)
1423 */
1424 for (i=0; i<mi->nmidas; i++) {
1425 mt = &mi->mterms[i];
1426 if (beta_type(mt->type)) {
1427 bound_rows += 2;
1428 } else if (mt->type == MIDAS_NEALMON) {
1429 bound_rows += mt->nparm;
1430 }
1431 if (mt->type == MIDAS_U) {
1432 nb += mt->nlags;
1433 } else if (takes_coeff(mt->type)) {
1434 nb += 1;
1435 }
1436 }
1437
1438 mi->bounds = gretl_matrix_alloc(bound_rows, 3);
1439 if (mi->bounds == NULL) {
1440 return E_ALLOC;
1441 }
1442
1443 /* build X matrix for cond. OLS */
1444 mi->X = gretl_matrix_alloc(mi->nobs, nb);
1445 if (mi->X == NULL) {
1446 return E_ALLOC;
1447 }
1448
1449 targ = mi->X->val;
1450 for (i=0; i<mi->nx; i++) {
1451 vi = mi->list[i+2];
1452 src = dset->Z[vi] + dset->t1;
1453 memcpy(targ, src, mi->colsize);
1454 targ += mi->nobs;
1455 }
1456
1457 k = j = 0;
1458
1459 for (i=0; i<mi->nmidas; i++) {
1460 mt = &mi->mterms[i];
1461 if (mi->bounds != NULL) {
1462 if (beta_type(mt->type)) {
1463 /* set up beta minima and maxima */
1464 for (ii=0; ii<2; ii++) {
1465 /* columns: index, minimum, maximum */
1466 gretl_matrix_set(mi->bounds, j, 0, k + ii + 1);
1467 if (mt->type == MIDAS_BETA1) {
1468 if (ii == 0) {
1469 gretl_matrix_set(mi->bounds, j, 1, mt->theta->val[0]);
1470 gretl_matrix_set(mi->bounds, j, 2, mt->theta->val[0]);
1471 } else {
1472 gretl_matrix_set(mi->bounds, j, 1, 1.0 + eps);
1473 gretl_matrix_set(mi->bounds, j, 2, 30);
1474 }
1475 } else {
1476 gretl_matrix_set(mi->bounds, j, 1, eps);
1477 gretl_matrix_set(mi->bounds, j, 2, 500);
1478 }
1479 j++;
1480 }
1481 } else if (mt->type == MIDAS_NEALMON) {
1482 /* set up nealmon minima and maxima */
1483 for (ii=0; ii<mt->nparm; ii++) {
1484 gretl_matrix_set(mi->bounds, j, 0, k + ii + 1);
1485 gretl_matrix_set(mi->bounds, j, 1, -2.0);
1486 gretl_matrix_set(mi->bounds, j, 2, +2.0);
1487 j++;
1488 }
1489 }
1490 }
1491 if (mt->type == MIDAS_U) {
1492 /* transcribe HF lags data */
1493 for (ii=0; ii<mt->nlags; ii++) {
1494 vi = mt->laglist[ii+1];
1495 src = dset->Z[vi] + dset->t1;
1496 memcpy(targ, src, mi->colsize);
1497 targ += mi->nobs;
1498 }
1499 } else {
1500 k += mt->nparm;
1501 if (takes_coeff(mt->type)) {
1502 /* leave space for weighted combo */
1503 targ += mi->nobs;
1504 }
1505 }
1506 }
1507
1508 mi->g = gretl_column_vector_alloc(k);
1509 mi->b = gretl_column_vector_alloc(nb);
1510 mi->theta = gretl_column_vector_alloc(k);
1511
1512 if (mi->g == NULL || mi->b == NULL || mi->theta == NULL) {
1513 return E_ALLOC;
1514 }
1515
1516 type0 = mi->mterms[0].type;
1517
1518 if (mi->nmidas == 1 && (type0 == MIDAS_BETA0 || type0 == MIDAS_BETA1)) {
1519 /* FIXME multiple beta terms? */
1520 midas_beta_init(mi);
1521 } else {
1522 /* initialize overall mi->theta by composition
1523 of one or more individual theta vectors
1524 */
1525 ii = 0;
1526 for (i=0; i<mi->nmidas; i++) {
1527 mt = &mi->mterms[i];
1528 if (mt->type != MIDAS_U) {
1529 for (j=0; j<mt->nparm; j++) {
1530 mi->theta->val[ii++] = mt->theta->val[j];
1531 }
1532 }
1533 }
1534 }
1535
1536 return 0;
1537 }
1538
midas_info_destroy(midas_info * mi)1539 static void midas_info_destroy (midas_info *mi)
1540 {
1541 gretl_matrix_free(mi->b);
1542 gretl_matrix_free(mi->g);
1543 gretl_matrix_free(mi->theta);
1544 gretl_matrix_free(mi->bounds);
1545 gretl_matrix_free(mi->u);
1546 gretl_matrix_free(mi->X);
1547 gretl_matrix_free(mi->y);
1548 gretl_matrix_free(mi->xi);
1549
1550 free(mi->mterms);
1551 free(mi->xlist);
1552 free(mi->seplist);
1553 g_free(mi->pnames);
1554
1555 free(mi);
1556 }
1557
1558 /* combined callback that returns SSR given @theta, and also
1559 calculates the gradient in @g
1560 */
1561
cond_ols_callback(double * theta,double * g,int n,void * ptr)1562 static double cond_ols_callback (double *theta, double *g,
1563 int n, void *ptr)
1564 {
1565 midas_info *mi = ptr;
1566 DATASET *dset = mi->dset;
1567 gretl_matrix *w;
1568 gretl_matrix *mg;
1569 int i, j, k, t, s;
1570 int vj, xcol;
1571 int err = 0;
1572
1573 #if MIDAS_DEBUG > 1
1574 fprintf(stderr, "cond_ols_callback: starting\n");
1575 for (i=0; i<n; i++) {
1576 fprintf(stderr, " theta[%d] = %g\n", i, theta[i]);
1577 }
1578 #endif
1579
1580 if (mi->nalmonp > 0) {
1581 /* reset mi->y, it may have been altered */
1582 memcpy(mi->y->val, dset->Z[mi->yno] + dset->t1, mi->colsize);
1583 }
1584
1585 /* Fill relevant columns of mi->X with weighted
1586 linear combinations of HF data
1587 */
1588
1589 /* initialize index of column to fill */
1590 xcol = mi->nx;
1591
1592 /* initialize index for reading from @theta */
1593 k = 0;
1594
1595 for (i=0; i<mi->nmidas && !err; i++) {
1596 midas_term *mt = &mi->mterms[i];
1597
1598 if (mt->type == MIDAS_U) {
1599 /* copy X data: job done already */
1600 xcol += mt->nlags;
1601 continue;
1602 }
1603 for (j=0; j<mt->nparm; j++) {
1604 mt->theta->val[j] = theta[k++];
1605 }
1606 w = midas_weights(mt->nlags, mt->theta, mt->type, &err);
1607 if (!err && mt->type == MIDAS_ALMONP) {
1608 /* net the estimated effect out of @y */
1609 for (j=0; j<mt->nlags; j++) {
1610 vj = mt->laglist[j+1];
1611 s = 0;
1612 for (t=dset->t1; t<=dset->t2; t++) {
1613 mi->y->val[s++] -= dset->Z[vj][t] * w->val[j];
1614 }
1615 }
1616 } else if (!err) {
1617 /* fill the next X column and advance */
1618 gretl_matrix_zero(mi->xi);
1619 for (j=0; j<mt->nlags; j++) {
1620 vj = mt->laglist[j+1];
1621 s = 0;
1622 for (t=dset->t1; t<=dset->t2; t++) {
1623 mi->xi->val[s++] += dset->Z[vj][t] * w->val[j];
1624 }
1625 }
1626 for (t=0; t<mi->nobs; t++) {
1627 gretl_matrix_set(mi->X, t, xcol, mi->xi->val[t]);
1628 }
1629 xcol++;
1630 }
1631 gretl_matrix_free(w);
1632 }
1633
1634 if (!err) {
1635 /* run OLS conditional on theta */
1636 err = gretl_matrix_ols(mi->y, mi->X, mi->b, NULL, mi->u, NULL);
1637 }
1638
1639 if (g == NULL) {
1640 goto skip_gradient;
1641 }
1642
1643 /* now for the gradient with respect to the hyper-parameters
1644 in @theta
1645 */
1646
1647 k = 0;
1648 xcol = mi->nx;
1649
1650 for (i=0; i<mi->nmidas && !err; i++) {
1651 midas_term *mt = &mi->mterms[i];
1652 double xit;
1653 int p, vp;
1654
1655 if (mt->type == MIDAS_U) {
1656 xcol += mt->nlags;
1657 continue;
1658 }
1659 mg = midas_gradient(mt->nlags, mt->theta, mt->type, &err);
1660 if (err) {
1661 break;
1662 }
1663 for (j=0; j<mg->cols; j++) {
1664 /* loop across hyper-parameters for this term */
1665 gretl_matrix_zero(mi->xi);
1666 for (t=0; t<mi->nobs; t++) {
1667 /* loop across observations */
1668 s = t + dset->t1;
1669 xit = 0.0;
1670 for (p=0; p<mt->nlags; p++) {
1671 /* loop across HF lags at this observation */
1672 vp = mt->laglist[p+1];
1673 xit += dset->Z[vp][s] * gretl_matrix_get(mg, p, j);
1674 }
1675 mi->xi->val[t] = xit;
1676 }
1677 if (takes_coeff(mt->type)) {
1678 gretl_matrix_multiply_by_scalar(mi->xi, mi->b->val[xcol]);
1679 }
1680 /* compute gradient of SSR wrt hyper-parameter @k */
1681 g[k] = 0.0;
1682 for (t=0; t<mi->nobs; t++) {
1683 g[k] -= 2 * mi->xi->val[t] * mi->u->val[t];
1684 }
1685 /* advance to next hyper-parameter */
1686 k++;
1687 }
1688 if (takes_coeff(mt->type)) {
1689 /* advance read position for OLS coeff */
1690 xcol++;
1691 }
1692 gretl_matrix_free(mg);
1693 }
1694
1695 skip_gradient:
1696
1697 if (err) {
1698 const char *s = errmsg_get_with_default(err);
1699
1700 fprintf(stderr, "cond_ols_callback: err=%d (%s)\n", err, s);
1701 return NADBL;
1702 }
1703
1704 mi->SSR = 0.0;
1705 for (t=0; t<mi->nobs; t++) {
1706 mi->SSR += mi->u->val[t] * mi->u->val[t];
1707 }
1708
1709 return mi->SSR;
1710 }
1711
1712 /* Given mi->b (OLS coefficients) and mi->theta (hyper-
1713 parameters optimized by L-BFGS-B), combine them into
1714 the full coefficient array.
1715 */
1716
make_full_coeff_vector(midas_info * mi,int n)1717 static int make_full_coeff_vector (midas_info *mi, int n)
1718 {
1719 gretl_matrix *c;
1720 midas_term *mt;
1721 int i, j, k, l, m;
1722
1723 c = gretl_column_vector_alloc(n);
1724
1725 if (c == NULL) {
1726 return E_ALLOC;
1727 }
1728
1729 k = l = m = 0;
1730
1731 for (i=0; i<mi->nx; i++) {
1732 c->val[k++] = mi->b->val[l++];
1733 }
1734 for (i=0; i<mi->nmidas; i++) {
1735 mt = &mi->mterms[i];
1736 if (mt->type == MIDAS_U) {
1737 for (j=0; j<mt->nlags; j++) {
1738 c->val[k++] = mi->b->val[l++];
1739 }
1740 } else {
1741 if (takes_coeff(mt->type)) {
1742 c->val[k++] = mi->b->val[l++];
1743 }
1744 for (j=0; j<mt->nparm; j++) {
1745 c->val[k++] = mi->theta->val[m++];
1746 }
1747 }
1748 }
1749
1750 /* swap out the original mi->b */
1751 gretl_matrix_free(mi->b);
1752 mi->b = c;
1753
1754 return 0;
1755 }
1756
transcribe_to_nlspec(nlspec * s,midas_info * mi,gretlopt opt)1757 static int transcribe_to_nlspec (nlspec *s,
1758 midas_info *mi,
1759 gretlopt opt)
1760 {
1761 int err = 0;
1762
1763 s->ncoeff = mi->b->rows + mi->theta->rows;
1764 err = make_full_coeff_vector(mi, s->ncoeff);
1765 if (err) {
1766 return err;
1767 }
1768
1769 s->ci = MIDASREG;
1770 s->opt = opt;
1771 s->dv = mi->list[1];
1772 s->coeff = mi->b->val;
1773 s->t1 = mi->dset->t1;
1774 s->t2 = mi->dset->t2;
1775 s->crit = mi->SSR;
1776 s->parnames = mi->pnames;
1777 s->fvec = mi->u->val;
1778 s->dset = mi->dset;
1779
1780 return err;
1781 }
1782
1783 /* The Gauss-Newton Regression after running BFGS plus
1784 conditional OLS */
1785
cond_ols_GNR(MODEL * pmod,midas_info * mi,gretlopt opt,PRN * prn)1786 static int cond_ols_GNR (MODEL *pmod,
1787 midas_info *mi,
1788 gretlopt opt,
1789 PRN *prn)
1790 {
1791 DATASET *dset = mi->dset;
1792 DATASET *gdset = NULL;
1793 int *glist = NULL;
1794 gretl_matrix *w, *G;
1795 int nc, zcol, vi;
1796 int i, j, k, kk, t, s, v;
1797 int cpi, *cpos = NULL;
1798 int err = 0;
1799
1800 #if MIDAS_DEBUG
1801 fprintf(stderr, "cond_ols_GNR, starting\n");
1802 #endif
1803
1804 /* total number of coefficients in the GNR */
1805 nc = mi->b->rows + mi->theta->rows;
1806
1807 v = nc + 2;
1808 for (i=2; i<=mi->list[0]; i++) {
1809 if (mi->list[i] == 0) {
1810 v--;
1811 break;
1812 }
1813 }
1814
1815 glist = gretl_list_new(nc + 1);
1816
1817 gdset = create_auxiliary_dataset(v, mi->nobs, 0);
1818 if (gdset == NULL) {
1819 return E_ALLOC;
1820 }
1821
1822 if (dataset_is_time_series(dset)) {
1823 gdset->structure = dset->structure;
1824 gdset->pd = dset->pd;
1825 ntolabel(gdset->stobs, dset->t1, dset);
1826 gdset->sd0 = get_date_x(gdset->pd, gdset->stobs);
1827 }
1828
1829 /* dependent var: nls residual */
1830 glist[0] = glist[1] = 1;
1831 memcpy(gdset->Z[1], mi->u->val, mi->colsize);
1832 strcpy(gdset->varname[1], "gy");
1833
1834 zcol = 2;
1835
1836 /* low-frequency regressors */
1837 for (i=2; i<=mi->list[0]; i++) {
1838 int gvi;
1839
1840 vi = mi->list[i];
1841 if (vi == 0) {
1842 gvi = 0;
1843 } else {
1844 gvi = zcol++;
1845 memcpy(gdset->Z[gvi], dset->Z[vi] + dset->t1, mi->colsize);
1846 sprintf(gdset->varname[gvi], "gx%d", i-1);
1847 }
1848 glist[i] = gvi;
1849 glist[0] += 1;
1850 }
1851
1852 if (mi->nbeta1 > 0) {
1853 cpos = gretl_list_new(mi->nbeta1);
1854 if (cpos == NULL) {
1855 err = E_ALLOC;
1856 }
1857 }
1858
1859 /* the read position for OLS coefficients */
1860 k = mi->nx;
1861
1862 /* the write position and coeff position for beta1 terms */
1863 cpi = 1;
1864 kk = mi->nx;
1865
1866 /* MIDAS terms */
1867 for (i=0; i<mi->nmidas && !err; i++) {
1868 midas_term *mt = &mi->mterms[i];
1869 double zt, gij, hfb = 1.0;
1870 int ii, pos;
1871
1872 if (mt->type == MIDAS_U) {
1873 /* gradient wrt U-MIDAS coeffs */
1874 for (j=0; j<mt->nlags; j++) {
1875 vi = mt->laglist[j+1];
1876 memcpy(gdset->Z[zcol], dset->Z[vi] + dset->t1,
1877 mi->colsize);
1878 glist[0] += 1;
1879 glist[glist[0]] = zcol;
1880 sprintf(gdset->varname[zcol], "mdx%d", i+1);
1881 zcol++;
1882 }
1883 k += mt->nlags;
1884 kk += mt->nlags;
1885 continue;
1886 }
1887
1888 /* gradient wrt weighted linear combination */
1889 w = midas_weights(mt->nlags, mt->theta, mt->type, &err);
1890 if (!err) {
1891 s = 0;
1892 for (t=dset->t1; t<=dset->t2; t++) {
1893 zt = 0.0;
1894 for (j=0; j<mt->nlags; j++) {
1895 vi = mt->laglist[j+1];
1896 zt += dset->Z[vi][t] * w->val[j];
1897 }
1898 gdset->Z[zcol][s] = zt;
1899 s++;
1900 }
1901 gretl_matrix_free(w);
1902 }
1903 glist[0] += 1;
1904 glist[glist[0]] = zcol;
1905 sprintf(gdset->varname[zcol], "mdx%d", i+1);
1906 zcol++;
1907 kk++;
1908
1909 /* gradient wrt hyper-parameters */
1910 pos = zcol;
1911 if (takes_coeff(mt->type)) {
1912 hfb = mi->b->val[k++];
1913 }
1914 G = midas_gradient(mt->nlags, mt->theta, mt->type, &err);
1915 if (!err) {
1916 s = 0;
1917 for (t=dset->t1; t<=dset->t2; t++) {
1918 for (ii=0; ii<mt->nparm; ii++) {
1919 zt = 0.0;
1920 for (j=0; j<mt->nlags; j++) {
1921 vi = mt->laglist[j+1];
1922 gij = gretl_matrix_get(G, j, ii);
1923 zt += dset->Z[vi][t] * gij;
1924 }
1925 gdset->Z[pos+ii][s] = zt * hfb;
1926 }
1927 s++;
1928 }
1929 gretl_matrix_free(G);
1930 }
1931
1932 if (!err) {
1933 for (ii=0; ii<mt->nparm; ii++) {
1934 glist[0] += 1;
1935 glist[glist[0]] = pos + ii;
1936 if (ii == 0 && mt->type == MIDAS_BETA1) {
1937 cpos[cpi++] = kk;
1938 sprintf(gdset->varname[pos+ii], "grad%dc", ii+1);
1939 } else {
1940 sprintf(gdset->varname[pos+ii], "grad%d", ii+1);
1941 }
1942 }
1943 zcol += mt->nparm;
1944 kk += mt->nparm;
1945 }
1946 }
1947
1948 #if MIDAS_DEBUG
1949 printlist(glist, "glist, mi");
1950 printlist(cpos, "cpos");
1951 # if MIDAS_DEBUG > 1
1952 printdata(glist, NULL, gdset, OPT_O, prn);
1953 # endif
1954 #endif
1955
1956 *pmod = GNR(glist, gdset, opt, prn);
1957
1958 #if MIDAS_DEBUG
1959 printmodel(pmod, gdset, opt, prn);
1960 #endif
1961
1962 if (pmod->errcode == 0 || pmod->errcode == E_JACOBIAN) {
1963 nlspec spec = {0};
1964
1965 if (pmod->errcode == 0 && (opt & OPT_B)) {
1966 /* test for structural break */
1967 QLR_test(pmod, gdset, OPT_Q | OPT_M, NULL);
1968 }
1969 err = transcribe_to_nlspec(&spec, mi, opt);
1970 if (!err) {
1971 err = finalize_nls_model(pmod, &spec, 0, glist);
1972 if (!err && cpos != NULL) {
1973 /* invalidate std errors of clamped coeffs */
1974 for (cpi=1; cpi<=cpos[0]; cpi++) {
1975 i = cpos[cpi];
1976 if (i < pmod->ncoeff) {
1977 /* this conditionality should NOT be needed! */
1978 pmod->sderr[i] = NADBL;
1979 }
1980 }
1981 }
1982 }
1983 if (err && !pmod->errcode) {
1984 pmod->errcode = err;
1985 }
1986 if (!pmod->errcode) {
1987 set_model_id(pmod, opt);
1988 }
1989 }
1990
1991 destroy_dataset(gdset);
1992 free(glist);
1993 free(cpos);
1994
1995 return err;
1996 }
1997
1998 /* Golden Section search for one-dimensional optimization:
1999 faster (though somewhat less robust?) than using L_BFGS-B.
2000 */
2001
midas_gss(double * theta,int n,midas_info * mi,int * fncount)2002 static int midas_gss (double *theta, int n, midas_info *mi,
2003 int *fncount)
2004 {
2005 double gr = (sqrt(5.0) + 1) / 2.0;
2006 double a = gretl_matrix_get(mi->bounds, 1, 1);
2007 double b = gretl_matrix_get(mi->bounds, 1, 2);
2008 double c = b - (b - a) / gr;
2009 double d = a + (b - a) / gr;
2010 double tol = 1.0e-4;
2011 double fc, fd;
2012 int err = 0;
2013
2014 #if MIDAS_DEBUG
2015 fprintf(stderr, "gss: bounds = [%g,%g], incoming value %g\n",
2016 a, b, theta[1]);
2017 #endif
2018
2019 /* theta[0] will be the clamped term, and
2020 theta[1] the one we're optimizing
2021 */
2022
2023 while (fabs(c - d) > tol) {
2024 theta[1] = c;
2025 fc = cond_ols_callback(theta, NULL, n, mi);
2026 theta[1] = d;
2027 fd = cond_ols_callback(theta, NULL, n, mi);
2028 *fncount += 2;
2029 if (na(fc) || na(fd)) {
2030 err = E_NAN;
2031 break;
2032 }
2033 if (fc < fd) {
2034 b = d;
2035 } else {
2036 a = c;
2037 }
2038 /* recompute c and d to preserve precision */
2039 c = b - (b - a) / gr;
2040 d = a + (b - a) / gr;
2041 }
2042
2043 if (!err) {
2044 /* record the optimum */
2045 theta[1] = (b + a) / 2;
2046 }
2047
2048 return err;
2049 }
2050
2051 /* driver function for estimating MIDAS model using
2052 L-BFGS-B or Golden Section search */
2053
midas_bfgs_run(MODEL * pmod,midas_info * mi,gretlopt opt,PRN * prn)2054 static int midas_bfgs_run (MODEL *pmod, midas_info *mi,
2055 gretlopt opt, PRN *prn)
2056 {
2057 int n = gretl_vector_get_length(mi->theta);
2058 double *theta = mi->theta->val;
2059 int fncount = 0, grcount = 0;
2060 int err;
2061
2062 if (theta[0] == 1.0 && mi->nmidas == 1 && mi->nbeta1 == 1) {
2063 #if MIDAS_DEBUG
2064 fprintf(stderr, "midas_bfgs_run: calling midas_gss\n");
2065 #endif
2066 err = midas_gss(theta, n, mi, &fncount);
2067 } else {
2068 double reltol = libset_get_double(BFGS_TOLER);
2069
2070 #if MIDAS_DEBUG
2071 fprintf(stderr, "midas_bfgs_run: calling LBFGS_max\n");
2072 #endif
2073 err = LBFGS_max(theta, n, 1000, reltol,
2074 &fncount, &grcount, NULL,
2075 C_SSR, NULL, cond_ols_callback,
2076 mi, mi->bounds, opt, prn);
2077 }
2078
2079 if (!err) {
2080 err = cond_ols_GNR(pmod, mi, opt, prn);
2081 }
2082
2083 if (!err) {
2084 gretl_model_set_int(pmod, "iters", fncount);
2085 }
2086
2087 return err;
2088 }
2089
2090 /* Below: we add a column vector containing the "gross"
2091 MIDAS coefficients (i.e. high-frequency slope, where
2092 applicable, times the weights implied by the hyper-
2093 parameter values). This is required for forecasting.
2094 In addition, if there's just a single MIDAS term, we
2095 add to the model a two-column matrix holding the
2096 gross coefficients and the associated lags -- this is
2097 to support a plot of the coefficients in the GUI.
2098 */
2099
add_gross_midas_coeffs(MODEL * pmod,midas_info * mi)2100 static int add_gross_midas_coeffs (MODEL *pmod,
2101 midas_info *mi)
2102 {
2103 const double *b = pmod->coeff + mi->nx;
2104 gretl_matrix *hfb, *C = NULL;
2105 int j, k = 0, nt = 0;
2106 int row = 0;
2107 int err = 0;
2108
2109 for (j=0; j<mi->nmidas; j++) {
2110 nt += mi->mterms[j].nlags;
2111 }
2112
2113 /* vector to hold all gross midas coeffs */
2114 hfb = gretl_column_vector_alloc(nt);
2115 if (hfb == NULL) {
2116 err = E_ALLOC;
2117 }
2118
2119 if (!err && mi->nmidas == 1) {
2120 /* matrix for plotting purposes */
2121 C = gretl_matrix_alloc(mi->mterms[0].nlags, 2);
2122 }
2123
2124 for (j=0; j<mi->nmidas && !err; j++) {
2125 midas_term *mt = &mi->mterms[j];
2126 gretl_matrix *theta = NULL;
2127
2128 if (mt->type != MIDAS_U) {
2129 theta = gretl_matrix_alloc(mt->nparm, 1);
2130 if (theta == NULL) {
2131 err = E_ALLOC;
2132 }
2133 }
2134
2135 if (!err) {
2136 gretl_matrix *w = NULL;
2137 double ci, hfslope = 0;
2138 int p = mt->minlag;
2139 int i;
2140
2141 if (mt->type == MIDAS_U) {
2142 for (i=0; i<mt->nparm; i++) {
2143 gretl_vector_set(hfb, row++, b[k]);
2144 if (C != NULL) {
2145 gretl_matrix_set(C, i, 0, b[k]);
2146 gretl_matrix_set(C, i, 1, p++);
2147 }
2148 k++;
2149 }
2150 } else {
2151 /* compute slope times weights */
2152 hfslope = takes_coeff(mt->type) ? b[k++] : 1.0;
2153 for (i=0; i<mt->nparm; i++) {
2154 theta->val[i] = b[k++];
2155 }
2156 w = midas_weights(mt->nlags, theta, mt->type, &err);
2157 if (!err) {
2158 for (i=0; i<mt->nlags; i++) {
2159 ci = hfslope * w->val[i];
2160 gretl_vector_set(hfb, row++, ci);
2161 if (C != NULL) {
2162 gretl_matrix_set(C, i, 0, ci);
2163 gretl_matrix_set(C, i, 1, p++);
2164 }
2165 }
2166 }
2167 gretl_matrix_free(w);
2168 }
2169 }
2170
2171 if (theta != NULL) {
2172 gretl_matrix_free(theta);
2173 theta = NULL;
2174 }
2175
2176 if (!err && C != NULL) {
2177 /* finalize plotting matrix and attach to model */
2178 char *cnames[2] = {"coeff", "lag"};
2179 char **S = strings_array_dup(cnames, 2);
2180
2181 if (S != NULL) {
2182 gretl_matrix_set_colnames(C, S);
2183 }
2184 gretl_model_set_matrix_as_data(pmod, "midas_coeffs", C);
2185 } else {
2186 gretl_matrix_free(C);
2187 }
2188 }
2189
2190 if (!err) {
2191 gretl_model_set_matrix_as_data(pmod, "hfb", hfb);
2192 }
2193
2194 return err;
2195 }
2196
2197 /* Add to @pmod an array of bundles holding information
2198 on the individual MIDAS terms (although in most cases
2199 there will be just one such bundle).
2200 */
2201
model_add_mterms_array(MODEL * pmod,midas_term * mterms,int nmidas)2202 static int model_add_mterms_array (MODEL *pmod,
2203 midas_term *mterms,
2204 int nmidas)
2205 {
2206 gretl_array *A;
2207 int err = 0;
2208
2209 A = gretl_array_new(GRETL_TYPE_BUNDLES, nmidas, &err);
2210
2211 if (A != NULL) {
2212 midas_term *mt;
2213 gretl_bundle *b;
2214 int i;
2215
2216 for (i=0; i<nmidas && !err; i++) {
2217 b = gretl_bundle_new();
2218 if (b == NULL) {
2219 err = E_ALLOC;
2220 } else {
2221 mt = &mterms[i];
2222 err += gretl_bundle_set_string(b, "lname", mt->lnam0);
2223 err += gretl_bundle_set_int(b, "prelag", prelag(mt) ? 1 : 0);
2224 err += gretl_bundle_set_int(b, "minlag", mt->minlag);
2225 err += gretl_bundle_set_int(b, "maxlag", mt->maxlag);
2226 err += gretl_bundle_set_int(b, "type", mt->type);
2227 err += gretl_bundle_set_int(b, "nparm", mt->nparm);
2228 err += gretl_array_set_bundle(A, i, b, 0);
2229 if (err) {
2230 err = E_DATA;
2231 }
2232 }
2233 }
2234
2235 if (err) {
2236 gretl_array_destroy(A);
2237 } else {
2238 gretl_model_set_array_as_data(pmod, "midas_info", A);
2239 }
2240 }
2241
2242 return err;
2243 }
2244
2245 /* Retrieve from array @A, which was earlier attached to
2246 a model, information on the individual MIDAS terms
2247 for the purpose of forecasting.
2248 */
2249
mterms_from_array(gretl_array * A,int * nmidas,int * err)2250 static midas_term *mterms_from_array (gretl_array *A,
2251 int *nmidas,
2252 int *err)
2253 {
2254 int n = gretl_array_get_length(A);
2255 midas_term *mt, *mterms = NULL;
2256
2257 if (n == 0) {
2258 *err = E_DATA;
2259 } else {
2260 mterms = malloc(n * sizeof *mterms);
2261 if (mterms == NULL) {
2262 *err = E_ALLOC;
2263 }
2264 }
2265
2266 if (mterms != NULL) {
2267 gretl_bundle *b;
2268 int i;
2269
2270 for (i=0; i<n && !*err; i++) {
2271 mt = &mterms[i];
2272 midas_term_init(mt);
2273 b = gretl_array_get_bundle(A, i);
2274 if (b == NULL) {
2275 *err = E_DATA;
2276 } else {
2277 strcpy(mt->lname, gretl_bundle_get_string(b, "lname", err));
2278 if (gretl_bundle_get_int(b, "prelag", err)) {
2279 mt->flags |= M_PRELAG;
2280 }
2281 mt->minlag = gretl_bundle_get_int(b, "minlag", err);
2282 mt->maxlag = gretl_bundle_get_int(b, "maxlag", err);
2283 mt->type = gretl_bundle_get_int(b, "type", err);
2284 mt->nparm = gretl_bundle_get_int(b, "nparm", err);
2285 sprintf(mt->mname, "theta___%d", i+1);
2286 }
2287 }
2288
2289 if (*err) {
2290 free(mterms);
2291 mterms = NULL;
2292 }
2293 }
2294
2295 if (!*err) {
2296 *nmidas = n;
2297 }
2298
2299 return mterms;
2300 }
2301
2302 /* For use when printing a midasreg model: returns a line
2303 to be printed before the coefficients associated with a
2304 MIDAS term, or NULL on failure.
2305 */
2306
get_midas_term_line(const MODEL * pmod,int i)2307 char *get_midas_term_line (const MODEL *pmod, int i)
2308 {
2309 gretl_array *A = gretl_model_get_data(pmod, "midas_info");
2310 gretl_bundle *b = NULL;
2311 char *ret = NULL;
2312
2313 if (A != NULL) {
2314 int n = gretl_array_get_length(A);
2315
2316 if (i >= 0 && i < n) {
2317 b = gretl_array_get_bundle(A, i);
2318 }
2319 }
2320
2321 if (b != NULL) {
2322 const char *lname;
2323 int l1, l2, err = 0;
2324
2325 lname = gretl_bundle_get_string(b, "lname", &err);
2326 l1 = gretl_bundle_get_int(b, "minlag", &err);
2327 l2 = gretl_bundle_get_int(b, "maxlag", &err);
2328
2329 ret = gretl_strdup_printf(_("MIDAS list %s, high-frequency lags %d to %d"),
2330 lname, l1, l2);
2331 }
2332
2333 return ret;
2334 }
2335
2336 /* Get the MIDAS model ready for shipping out. What
2337 exactly we do here depends in part on whether
2338 estimation was done by NLS or OLS.
2339 */
2340
finalize_midas_model(MODEL * pmod,midas_info * mi,const char * param,gretlopt opt,const DATASET * dset)2341 static int finalize_midas_model (MODEL *pmod,
2342 midas_info *mi,
2343 const char *param,
2344 gretlopt opt,
2345 const DATASET *dset)
2346 {
2347 int type0, mixed = 0;
2348 int i, err = 0;
2349
2350 /* ensure that the incoming sample range is recorded */
2351 gretl_model_smpl_init(pmod, dset);
2352
2353 gretl_model_set_string_as_data(pmod, "midas_spec",
2354 gretl_strdup(param));
2355
2356 if (pmod->ci == OLS) {
2357 /* @pmod is the result of OLS estimation */
2358 int vi;
2359
2360 gretl_model_allocate_param_names(pmod, pmod->ncoeff);
2361 for (i=0; i<pmod->ncoeff; i++) {
2362 vi = pmod->list[i+2];
2363 gretl_model_set_param_name(pmod, i, dset->varname[vi]);
2364 }
2365 gretl_model_set_int(pmod, "umidas", 1);
2366 } else {
2367 /* @pmod is the result of some sort of NLS estimation */
2368 free(pmod->depvar);
2369 pmod->depvar = gretl_strdup(dset->varname[mi->list[1]]);
2370 free(pmod->list);
2371 pmod->list = gretl_list_copy(mi->list);
2372 if (mi->method == MDS_BFGS) {
2373 gretl_model_set_int(pmod, "BFGS", 1);
2374 } else if (mi->method == MDS_GSS) {
2375 gretl_model_set_int(pmod, "GSS", 1);
2376 }
2377 }
2378
2379 pmod->ci = MIDASREG;
2380
2381 /* record list of low-frequency regressors */
2382 if (mi->xlist == NULL) {
2383 gretl_model_set_int(pmod, "no_lfx", 1);
2384 } else {
2385 gretl_model_set_list_as_data(pmod, "lfxlist", mi->xlist);
2386 mi->xlist = NULL; /* donated to pmod */
2387 }
2388
2389 /* and positions where MIDAS terms start */
2390 if (mi->seplist != NULL) {
2391 gretl_model_set_list_as_data(pmod, "seplist", mi->seplist);
2392 mi->seplist = NULL; /* donated to pmod */
2393 }
2394
2395 if (mi->ldepvar) {
2396 gretl_model_set_int(pmod, "ldepvar", mi->ldepvar);
2397 gretl_model_set_int(pmod, "dynamic", 1);
2398 }
2399
2400 /* record the (common) type of MIDAS term? */
2401 type0 = mi->mterms[0].type;
2402 if (mi->nmidas > 1) {
2403 for (i=1; i<mi->nmidas; i++) {
2404 if (mi->mterms[i].type != type0) {
2405 mixed = 1;
2406 break;
2407 }
2408 }
2409 }
2410 if (!mixed && type0 > 0) {
2411 gretl_model_set_int(pmod, "midas_type", type0);
2412 }
2413
2414 if (!err) {
2415 /* Compute and record the "gross" MIDAS coefficients:
2416 we may need these for forecasting, and possibly to
2417 enable a plot in the GUI.
2418 */
2419 err = add_gross_midas_coeffs(pmod, mi);
2420 }
2421
2422 if (!err) {
2423 err = model_add_mterms_array(pmod, mi->mterms, mi->nmidas);
2424 }
2425
2426 return err;
2427 }
2428
2429 /* simple initialization of MIDAS slope terms for use with
2430 Levenberg-Marquardt */
2431
midas_means_init(midas_info * mi,const gretl_matrix * y,const gretl_matrix * X,gretl_matrix * c,const DATASET * dset)2432 static int midas_means_init (midas_info *mi,
2433 const gretl_matrix *y,
2434 const gretl_matrix *X,
2435 gretl_matrix *c,
2436 const DATASET *dset)
2437 {
2438 gretl_matrix *XZ;
2439 int err = 0;
2440
2441 XZ = gretl_matrix_alloc(mi->nobs, mi->nx + mi->hfslopes);
2442
2443 if (XZ == NULL) {
2444 err = E_ALLOC;
2445 } else {
2446 midas_term *mt;
2447 double zti;
2448 int i, j, k, s, t;
2449
2450 if (X != NULL) {
2451 /* transcribe low-frequency regressors */
2452 memcpy(XZ->val, X->val, mi->nobs * mi->nx * sizeof(double));
2453 }
2454
2455 /* transcribe time-mean of MIDAS terms */
2456 k = mi->nx * mi->nobs;
2457 for (i=0; i<mi->nmidas; i++) {
2458 mt = &mi->mterms[i];
2459 if (takes_coeff(mt->type)) {
2460 for (t=0; t<mi->nobs; t++) {
2461 s = t + dset->t1;
2462 zti = 0.0;
2463 for (j=0; j<mt->nlags; j++) {
2464 zti += dset->Z[mt->laglist[j+1]][s];
2465 }
2466 XZ->val[k++] = zti / mt->nlags;
2467 }
2468 }
2469 }
2470
2471 err = gretl_matrix_ols(y, XZ, c, NULL,
2472 NULL, NULL);
2473 gretl_matrix_free(XZ);
2474 }
2475
2476 return err;
2477 }
2478
2479 /* Define "private" matrices to hold the regular X data
2480 (MX___) and the vector of coefficients on these data
2481 (bx___). If we have an "hfslope" terms, then also try
2482 some initialization.
2483
2484 Note: we come here only if we're using native NLS
2485 (not U-MIDAS OLS, or BFGS + conditional OLS).
2486 */
2487
add_midas_matrices(midas_info * mi,const DATASET * dset)2488 static int add_midas_matrices (midas_info *mi,
2489 const DATASET *dset)
2490 {
2491 gretl_matrix *X = NULL;
2492 gretl_matrix *b = NULL;
2493 gretl_matrix *y = NULL;
2494 gretl_matrix *c = NULL;
2495 int init_err = 0;
2496 int err = 0;
2497
2498 if (mi->nx > 0) {
2499 X = gretl_matrix_data_subset(mi->xlist, dset,
2500 dset->t1, dset->t2,
2501 M_MISSING_ERROR,
2502 &err);
2503 if (!err) {
2504 err = private_matrix_add(X, "MX___");
2505 }
2506
2507 if (!err) {
2508 b = gretl_zero_matrix_new(mi->nx, 1);
2509 if (b != NULL) {
2510 err = private_matrix_add(b, "bx___");
2511 } else {
2512 err = E_ALLOC;
2513 }
2514 }
2515 }
2516
2517 if (mi->nx == 0 && mi->hfslopes == 0) {
2518 return err;
2519 }
2520
2521 if (!err) {
2522 /* for initialization only */
2523 y = gretl_column_vector_alloc(mi->nobs);
2524 if (y != NULL) {
2525 memcpy(y->val, dset->Z[mi->yno] + dset->t1,
2526 mi->colsize);
2527 } else {
2528 init_err = 1;
2529 }
2530 }
2531
2532 if (!err && !init_err) {
2533 /* "full-length" coeff vector */
2534 c = gretl_zero_matrix_new(mi->nx + mi->hfslopes, 1);
2535 if (c == NULL) {
2536 init_err = 1;
2537 }
2538 }
2539
2540 if (!err && !init_err) {
2541 if (mi->hfslopes > 0) {
2542 init_err = midas_means_init(mi, y, X, c, dset);
2543 }
2544 if (init_err || (mi->hfslopes == 0 && mi->nx > 0)) {
2545 /* simpler initialization, ignoring HF terms */
2546 c->rows = mi->nx;
2547 init_err = gretl_matrix_ols(y, X, c, NULL,
2548 NULL, NULL);
2549 }
2550 }
2551
2552 if (!err) {
2553 int i;
2554
2555 if (!init_err) {
2556 /* initialize X coeffs from OLS */
2557 for (i=0; i<mi->nx; i++) {
2558 b->val[i] = c->val[i];
2559 }
2560 }
2561 if (mi->hfslopes > 0) {
2562 /* initialize HF slopes, with fallback to zero */
2563 int use_c = !init_err && c->rows > mi->nx;
2564 char tmp[24];
2565 double bzi;
2566
2567 for (i=0; i<mi->nmidas && !err; i++) {
2568 if (takes_coeff(mi->mterms[i].type)) {
2569 sprintf(tmp, "bmlc___%d", i + 1);
2570 bzi = use_c ? c->val[mi->nx+i] : 0.0;
2571 err = private_scalar_add(bzi, tmp);
2572 }
2573 }
2574 }
2575 }
2576
2577 gretl_matrix_free(y);
2578 gretl_matrix_free(c);
2579
2580 return err;
2581 }
2582
put_midas_nls_line(char * line,DATASET * dset,gretlopt opt,PRN * prn)2583 static int put_midas_nls_line (char *line,
2584 DATASET *dset,
2585 gretlopt opt,
2586 PRN *prn)
2587 {
2588 int err;
2589
2590 if (opt & OPT_P) {
2591 /* display what we're passing to nls */
2592 pputs(prn, line);
2593 pputc(prn, '\n');
2594 }
2595
2596 err = nl_parse_line(NLS, line, dset, prn);
2597 *line = '\0';
2598
2599 return err;
2600 }
2601
2602 /* Append @pname to @s, ensuring that it's preceded by
2603 a single space if it's not preceded by a double quote.
2604 */
2605
append_pname(char * s,const char * pname)2606 static void append_pname (char *s, const char *pname)
2607 {
2608 if (*s != '\0') {
2609 char c = s[strlen(s) - 1];
2610
2611 if (c != ' ' && c != '"') {
2612 strcat(s, " ");
2613 }
2614 }
2615
2616 strcat(s, pname);
2617 }
2618
make_pname(char * targ,midas_term * mt,int i,const DATASET * dset)2619 static void make_pname (char *targ, midas_term *mt, int i,
2620 const DATASET *dset)
2621 {
2622 if (mt->type == MIDAS_NEALMON) {
2623 sprintf(targ, "Almon%d", i+1);
2624 } else if (beta_type(mt->type)) {
2625 sprintf(targ, "Beta%d", i+1);
2626 } else if (mt->type == MIDAS_ALMONP) {
2627 sprintf(targ, "Almon%d", i);
2628 } else {
2629 /* U-MIDAS */
2630 int *list = get_list_by_name(mt->lname);
2631
2632 if (list != NULL) {
2633 strcpy(targ, dset->varname[list[i+1]]);
2634 } else {
2635 sprintf(targ, "U-MIDAS%d", i+1);
2636 }
2637 }
2638 }
2639
add_param_names(midas_info * mi,const DATASET * dset)2640 static int add_param_names (midas_info *mi,
2641 const DATASET *dset)
2642 {
2643 char tmp[64], str[2*MAXLEN];
2644 int i, j, k = 0;
2645 int err = 0;
2646
2647 mi->seplist = gretl_list_new(mi->nmidas);
2648 *tmp = *str = '\0';
2649
2650 if (mi->xlist != NULL) {
2651 for (i=1; i<=mi->xlist[0]; i++) {
2652 strcpy(tmp, dset->varname[mi->xlist[i]]);
2653 append_pname(str, tmp);
2654 }
2655 k += mi->xlist[0];
2656 }
2657
2658 for (i=0; i<mi->nmidas; i++) {
2659 if (mi->seplist != NULL) {
2660 mi->seplist[i+1] = k;
2661 }
2662 if (takes_coeff(mi->mterms[i].type)) {
2663 if (mi->hfslopes > 1) {
2664 sprintf(tmp, "HF_slope%d", i+1);
2665 } else {
2666 strcpy(tmp, "HF_slope");
2667 }
2668 append_pname(str, tmp);
2669 k++;
2670 }
2671 for (j=0; j<mi->mterms[i].nparm; j++) {
2672 make_pname(tmp, &mi->mterms[i], j, dset);
2673 append_pname(str, tmp);
2674 k++;
2675 }
2676 }
2677
2678 mi->pnames = g_strdup(str);
2679
2680 return err;
2681 }
2682
2683 /* Allocate and initialize midas_info struct, and set
2684 MDS_NLS (Levenberg-Marquardt) as the default estimation
2685 method -- but this will be subject to revision when we
2686 digest the user's specification.
2687 */
2688
midas_info_new(const int * list,DATASET * dset)2689 static midas_info *midas_info_new (const int *list,
2690 DATASET *dset)
2691 {
2692 midas_info *mi = malloc(sizeof *mi);
2693
2694 if (mi != NULL) {
2695 /* zero everything first */
2696 memset(mi, 0, sizeof *mi);
2697 mi->yno = list[1];
2698 mi->list = list;
2699 mi->dset = dset;
2700 mi->method = MDS_NLS;
2701 mi->orig_toler = 0;
2702 }
2703
2704 return mi;
2705 }
2706
2707 /* When using Levenberg-Marquardt NLS, check to see if
2708 we have any terms that call for a "small" step in
2709 the optimizer (to try to head off numerical
2710 problems). This applies to the "highly nonlinear"
2711 parameterizations.
2712 */
2713
any_smallstep_terms(midas_info * mi)2714 static int any_smallstep_terms (midas_info *mi)
2715 {
2716 midas_term *mt;
2717 int i;
2718
2719 for (i=0; i<mi->nmidas; i++) {
2720 mt = &mi->mterms[i];
2721 if (mt->type == MIDAS_NEALMON ||
2722 mt->type == MIDAS_BETA0 ||
2723 mt->type == MIDAS_BETAN) {
2724 return 1;
2725 }
2726 }
2727
2728 return 0;
2729 }
2730
2731 /* General setup for NLS via Levenverg-Marquardt */
2732
midas_nls_setup(midas_info * mi,DATASET * dset,gretlopt opt,PRN * prn)2733 static int midas_nls_setup (midas_info *mi, DATASET *dset,
2734 gretlopt opt, PRN *prn)
2735 {
2736 char tmp[64], line[MAXLEN];
2737 midas_term *mt;
2738 const int *list = mi->list;
2739 int nmidas = mi->nmidas;
2740 int i, j = 0;
2741 int err = 0;
2742
2743 if (opt & OPT_P) {
2744 pputs(prn, "\n=== auto-generated nls specification ===\n");
2745 }
2746
2747 /* initial "nls" line */
2748 sprintf(line, "nls %s = ", dset->varname[list[1]]);
2749 if (mi->xlist != NULL) {
2750 strcat(line, "lincomb(XL___, bx___)");
2751 }
2752 for (i=0; i<nmidas; i++) {
2753 mt = &mi->mterms[i];
2754 if (takes_coeff(mt->type)) {
2755 sprintf(tmp, " + bmlc___%d*mlc___%d", ++j, i+1);
2756 } else {
2757 sprintf(tmp, " + mlc___%d", i+1);
2758 }
2759 strcat(line, tmp);
2760 }
2761 err = put_midas_nls_line(line, dset, opt, prn);
2762
2763 /* MIDAS series and gradient matrices */
2764 for (i=0; i<nmidas && !err; i++) {
2765 mt = &mi->mterms[i];
2766 if (mt->type == MIDAS_U) {
2767 sprintf(line, "series mlc___%d = lincomb(%s, %s)",
2768 i+1, mt->lname, mt->mname);
2769 } else {
2770 sprintf(line, "series mlc___%d = mlincomb(%s, %s, %d)",
2771 i+1, mt->lname, mt->mname, mt->type);
2772 }
2773 err = put_midas_nls_line(line, dset, opt, prn);
2774 if (!err && mt->type != MIDAS_U) {
2775 sprintf(line, "matrix mgr___%d = mgradient(%d, %s, %d)",
2776 i+1, mt->nlags, mt->mname, mt->type);
2777 err = put_midas_nls_line(line, dset, opt, prn);
2778 }
2779 }
2780
2781 /* derivatives */
2782 if (!err && mi->xlist != NULL) {
2783 strcpy(line, "deriv bx___ = MX___");
2784 err = put_midas_nls_line(line, dset, opt, prn);
2785 }
2786 for (i=0; i<nmidas && !err; i++) {
2787 mt = &mi->mterms[i];
2788 if (takes_coeff(mt->type)) {
2789 sprintf(line, "deriv bmlc___%d = {mlc___%d}", i+1, i+1);
2790 err = put_midas_nls_line(line, dset, opt, prn);
2791 }
2792 if (!err) {
2793 if (takes_coeff(mt->type)) {
2794 sprintf(line, "deriv %s = bmlc___%d * {%s} * mgr___%d",
2795 mt->mname, i+1, mt->lname, i+1);
2796 } else if (mt->type == MIDAS_ALMONP) {
2797 sprintf(line, "deriv %s = {%s} * mgr___%d",
2798 mt->mname, mt->lname, i+1);
2799 } else {
2800 sprintf(line, "deriv %s = {%s}", mt->mname, mt->lname);
2801 }
2802 err = put_midas_nls_line(line, dset, opt, prn);
2803 }
2804 }
2805
2806 if (!err) {
2807 /* add parameter names */
2808 sprintf(line, "param_names \"%s\"", mi->pnames);
2809 err = put_midas_nls_line(line, dset, opt, prn);
2810 }
2811
2812 if (opt & OPT_P) {
2813 pputs(prn, "=== end nls specification ===\n\n");
2814 }
2815
2816 if (!err) {
2817 if (any_smallstep_terms(mi)) {
2818 nl_set_smallstep();
2819 }
2820 if (mi->nmidas > 1) {
2821 /* allow for a sloppier than usual fit? */
2822 mi->orig_toler = libset_get_double(NLS_TOLER);
2823 libset_set_double(NLS_TOLER, 1.0e-7);
2824 }
2825 }
2826
2827 return err;
2828 }
2829
2830 /* Main driver function for built-in MIDAS estimation.
2831 The actual engine used is NLS, LBFGS or OLS.
2832 */
2833
midas_model(const int * list,const char * param,DATASET * dset,gretlopt opt,PRN * prn)2834 MODEL midas_model (const int *list,
2835 const char *param,
2836 DATASET *dset,
2837 gretlopt opt,
2838 PRN *prn)
2839 {
2840 midas_info *mi = NULL;
2841 int origv = dset->v;
2842 int save_t1 = dset->t1;
2843 int save_t2 = dset->t2;
2844 MODEL mod;
2845 int i, err = 0;
2846
2847 gretl_model_init(&mod, dset);
2848
2849 mi = midas_info_new(list, dset);
2850 if (mi == NULL) {
2851 mod.errcode = E_ALLOC;
2852 return mod;
2853 }
2854
2855 #if MIDAS_DEBUG
2856 fprintf(stderr, "\nmidas_model, starting...\n");
2857 #endif
2858
2859 if (param == NULL || *param == '\0') {
2860 err = E_DATA;
2861 } else {
2862 err = parse_midas_specs(mi, param, dset, opt);
2863 }
2864
2865 if (!err) {
2866 /* build list of regular regressors */
2867 err = make_midas_xlist(mi);
2868 if (mi->xlist != NULL &&
2869 (mi->method == MDS_NLS || mi->method == MDS_OLS)) {
2870 err = remember_list(mi->xlist, "XL___", NULL);
2871 user_var_privatize_by_name("XL___");
2872 }
2873 }
2874
2875 if (!err) {
2876 /* build (or just read) MIDAS lag-lists */
2877 err = make_midas_laglists(mi, dset);
2878 }
2879
2880 if (!err && mi->method == MDS_OLS) {
2881 err = umidas_ols(&mod, mi, dset, opt);
2882 goto midas_finish;
2883 }
2884
2885 if (!err) {
2886 /* determine usable sample range */
2887 err = midas_set_sample(mi, dset);
2888 }
2889
2890 if (!err && mi->method == MDS_NLS && mi->nx > 0) {
2891 /* add some required matrices */
2892 err = add_midas_matrices(mi, dset);
2893 }
2894
2895 if (!err) {
2896 /* construct string with names of parameters */
2897 err = add_param_names(mi, dset);
2898 }
2899
2900 if (!err && mi->method == MDS_NLS) {
2901 /* estimation using native NLS */
2902 err = midas_nls_setup(mi, dset, opt, prn);
2903 if (!err) {
2904 /* OPT_S: suppress gradient check */
2905 mod = nl_model(dset, (opt | OPT_S | OPT_M), prn);
2906 }
2907 } else if (!err) {
2908 /* estimation using L-BFGS-B or GSS */
2909 err = midas_bfgs_setup(mi, dset, opt);
2910 if (!err) {
2911 err = midas_bfgs_run(&mod, mi, opt, prn);
2912 }
2913 }
2914
2915 midas_finish:
2916
2917 dset->t1 = save_t1;
2918 dset->t2 = save_t2;
2919
2920 for (i=0; i<mi->nmidas; i++) {
2921 midas_term *mt = &mi->mterms[i];
2922
2923 if (!prelag(mt)) {
2924 free(mt->laglist);
2925 if (mt->flags & M_TMPLIST) {
2926 user_var_delete_by_name(mt->lname, NULL);
2927 }
2928 }
2929 if (mi->method == MDS_NLS && mt->type != MIDAS_U) {
2930 char tmp[24];
2931
2932 sprintf(tmp, "mgr___%d", i+1);
2933 user_var_delete_by_name(tmp, NULL);
2934 }
2935 }
2936
2937 if (err && !mod.errcode) {
2938 mod.errcode = err;
2939 }
2940
2941 if (!mod.errcode) {
2942 finalize_midas_model(&mod, mi, param, opt, dset);
2943 }
2944
2945 destroy_private_uvars();
2946
2947 if (mi != NULL) {
2948 if (mi->orig_toler != 0) {
2949 /* put back what we found on input */
2950 libset_set_double(NLS_TOLER, mi->orig_toler);
2951 }
2952 midas_info_destroy(mi);
2953 }
2954
2955 if (dset->v > origv) {
2956 /* or maybe not? */
2957 dataset_drop_last_variables(dset, dset->v - origv);
2958 }
2959
2960 return mod;
2961 }
2962