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 #include "matrix_extra.h"
23 #include "gretl_bfgs.h"
24 #include "libset.h"
25
26 #include <errno.h>
27
28 #define DDEBUG 0
29
30 typedef struct duration_info_ duration_info;
31
32 enum {
33 DUR_WEIBULL,
34 DUR_EXPON,
35 DUR_LOGLOG,
36 DUR_LOGNORM
37 };
38
39 enum {
40 DUR_UPDATE_XB = 1 << 0,
41 DUR_CONST_ONLY = 1 << 1
42 };
43
44 struct duration_info_ {
45 int dist; /* distribution type */
46 int flags; /* control info */
47 int k; /* number of covariates (including constant) */
48 int npar; /* total number of parameters */
49 int n; /* number of observations */
50 double ll; /* loglikelihood */
51 double *theta; /* parameter array, length npar */
52 gretl_matrix_block *B; /* workspace */
53 gretl_vector *logt; /* log of dependent variable (duration) */
54 gretl_matrix *X; /* covariates */
55 gretl_vector *cens; /* censoring variable (if needed) */
56 gretl_matrix *beta; /* coeffs on covariates */
57 gretl_matrix *llt; /* per-observation likelihood */
58 gretl_matrix *Xb; /* X \times \beta */
59 gretl_matrix *G; /* score */
60 gretl_matrix *V; /* covariance matrix */
61 PRN *prn; /* verbose printer */
62 };
63
duration_free(duration_info * dinfo)64 static void duration_free (duration_info *dinfo)
65 {
66 gretl_matrix_block_destroy(dinfo->B);
67 free(dinfo->theta);
68 }
69
70 /* initialize using OLS regression of the log of duration
71 on the covariates (or a simpler variant if we're
72 estimating the constant-only model)
73 */
74
duration_estimates_init(duration_info * dinfo)75 static int duration_estimates_init (duration_info *dinfo)
76 {
77 int err = 0;
78
79 if (dinfo->flags & DUR_CONST_ONLY) {
80 dinfo->theta[0] = gretl_vector_mean(dinfo->logt);
81 } else {
82 gretl_matrix *b = gretl_matrix_alloc(dinfo->k, 1);
83 int j;
84
85 if (b == NULL) {
86 return E_ALLOC;
87 }
88
89 err = gretl_matrix_ols(dinfo->logt, dinfo->X, b,
90 NULL, NULL, NULL);
91
92 if (!err) {
93 for (j=0; j<dinfo->k; j++) {
94 dinfo->theta[j] = b->val[j];
95 }
96 }
97
98 gretl_matrix_free(b);
99 }
100
101 if (dinfo->dist != DUR_EXPON) {
102 dinfo->theta[dinfo->k] = 1.0;
103 }
104
105 return err;
106 }
107
duration_init(duration_info * dinfo,MODEL * pmod,int censvar,const DATASET * dset,gretlopt opt,PRN * prn)108 static int duration_init (duration_info *dinfo, MODEL *pmod,
109 int censvar, const DATASET *dset,
110 gretlopt opt, PRN *prn)
111 {
112 int cn, n = pmod->nobs;
113 int np, k = pmod->ncoeff;
114 int i, j, t, v;
115 int err = 0;
116
117 dinfo->B = NULL;
118 dinfo->theta = NULL;
119
120 if (opt & OPT_E) {
121 /* exponential */
122 dinfo->dist = DUR_EXPON;
123 np = dinfo->npar = k;
124 } else if (opt & OPT_L) {
125 /* log-logistic */
126 dinfo->dist = DUR_LOGLOG;
127 np = dinfo->npar = k + 1;
128 } else if (opt & OPT_Z) {
129 /* log-normal */
130 dinfo->dist = DUR_LOGNORM;
131 np = dinfo->npar = k + 1;
132 } else {
133 /* default: Weibull */
134 dinfo->dist = DUR_WEIBULL;
135 np = dinfo->npar = k + 1;
136 }
137
138 dinfo->flags = 0;
139
140 dinfo->theta = malloc(np * sizeof *dinfo->theta);
141
142 if (dinfo->theta == NULL) {
143 return E_ALLOC;
144 }
145
146 cn = (censvar > 0)? n : 0;
147
148 dinfo->B = gretl_matrix_block_new(&dinfo->logt, n, 1,
149 &dinfo->X, n, k,
150 &dinfo->cens, cn, 1,
151 &dinfo->beta, k, 1,
152 &dinfo->Xb, n, 1,
153 &dinfo->llt, n, 1,
154 &dinfo->G, n, np,
155 NULL);
156 if (dinfo->B == NULL) {
157 return E_ALLOC;
158 }
159
160 if (cn == 0) {
161 /* mask unused zero-size part of matrix block */
162 dinfo->cens = NULL;
163 }
164
165 /* transcribe data into matrix form, taking the
166 log of the duration measurements */
167
168 i = 0;
169 for (t=pmod->t1; t<=pmod->t2; t++) {
170 if (na(pmod->uhat[t])) {
171 continue;
172 }
173 v = pmod->list[1];
174 dinfo->logt->val[i] = log(dset->Z[v][t]);
175 if (dinfo->cens != NULL) {
176 dinfo->cens->val[i] = dset->Z[censvar][t];
177 }
178 for (j=0; j<k; j++) {
179 v = pmod->list[j+2];
180 gretl_matrix_set(dinfo->X, i, j, dset->Z[v][t]);
181 }
182 i++;
183 }
184
185 dinfo->k = k;
186 dinfo->n = n;
187
188 err = duration_estimates_init(dinfo);
189
190 if (!err) {
191 dinfo->ll = NADBL;
192 dinfo->prn = (opt & OPT_V)? prn : NULL;
193 }
194
195 return err;
196 }
197
duration_update_Xb(duration_info * dinfo,const double * theta)198 static void duration_update_Xb (duration_info *dinfo, const double *theta)
199 {
200 int j;
201
202 if (theta == NULL) {
203 theta = dinfo->theta;
204 }
205
206 for (j=0; j<dinfo->k; j++) {
207 dinfo->beta->val[j] = theta[j];
208 }
209
210 gretl_matrix_multiply(dinfo->X, dinfo->beta, dinfo->Xb);
211 }
212
213 #define uncensored(d,i) (d->cens == NULL || d->cens->val[i] == 0)
214
215 /* The approach taken here to the loglikelihood and score for duration
216 models is that of Kalbfleisch and Prentice: see their Statistical
217 Analysis of Failure Time Data, 2e (Wiley, 2002), pp. 68-70.
218 */
219
duration_loglik(const double * theta,void * data)220 static double duration_loglik (const double *theta, void *data)
221 {
222 duration_info *dinfo = (duration_info *) data;
223 double *ll = dinfo->llt->val;
224 double *Xb = dinfo->Xb->val;
225 double *logt = dinfo->logt->val;
226 double wi, s = 1.0, lns = 0.0;
227 double l1ew = 0.0;
228 int i, di;
229
230 if (dinfo->dist != DUR_EXPON) {
231 s = theta[dinfo->k];
232 if (s <= 0) {
233 return NADBL;
234 }
235 lns = log(s);
236 }
237
238 duration_update_Xb(dinfo, theta);
239
240 dinfo->ll = 0.0;
241 errno = 0;
242
243 for (i=0; i<dinfo->n; i++) {
244 di = uncensored(dinfo, i);
245 wi = (logt[i] - Xb[i]) / s;
246 if (dinfo->dist == DUR_LOGLOG) {
247 l1ew = log(1 + exp(wi));
248 ll[i] = -l1ew;
249 if (di) {
250 ll[i] += wi - l1ew - lns;
251 }
252 } else if (dinfo->dist == DUR_LOGNORM) {
253 if (di) {
254 /* density */
255 ll[i] = -lns + log_normal_pdf(wi);
256 } else {
257 /* survivor */
258 ll[i] = log(normal_cdf(-wi));
259 }
260 } else {
261 /* Weibull, exponential */
262 ll[i] = -exp(wi);
263 if (di) {
264 ll[i] += wi - lns;
265 }
266 }
267 dinfo->ll += ll[i];
268 }
269
270 if (errno) {
271 dinfo->ll = NADBL;
272 }
273
274 return dinfo->ll;
275 }
276
277 /* normal hazard: ratio of density to survivor function */
278
normal_h(double w)279 static double normal_h (double w)
280 {
281 return normal_pdf(w) / normal_cdf(-w);
282 }
283
duration_score(double * theta,double * g,int np,BFGS_CRIT_FUNC ll,void * data)284 static int duration_score (double *theta, double *g, int np,
285 BFGS_CRIT_FUNC ll, void *data)
286 {
287 duration_info *dinfo = (duration_info *) data;
288 const double *logt = dinfo->logt->val;
289 const double *Xb = dinfo->Xb->val;
290 double wi, ewi, ai, xij, gij, s = 1.0;
291 int i, j, di, err = 0;
292
293 if (dinfo->flags == DUR_UPDATE_XB) {
294 duration_update_Xb(dinfo, theta);
295 }
296
297 if (dinfo->dist != DUR_EXPON) {
298 s = theta[dinfo->k];
299 }
300
301 if (g != NULL) {
302 for (j=0; j<np; j++) {
303 g[j] = 0.0;
304 }
305 }
306
307 for (i=0; i<dinfo->n; i++) {
308 di = uncensored(dinfo, i);
309 wi = (logt[i] - Xb[i]) / s;
310 ewi = exp(wi);
311 if (dinfo->dist == DUR_LOGLOG) {
312 ai = -di + (1 + di) * ewi / (1 + ewi);
313 } else if (dinfo->dist == DUR_LOGNORM) {
314 ai = di ? wi : normal_h(wi);
315 } else {
316 ai = ewi - di;
317 }
318 for (j=0; j<np; j++) {
319 if (j < dinfo->k) {
320 /* covariates */
321 xij = gretl_matrix_get(dinfo->X, i, j);
322 gij = xij * ai;
323 } else {
324 /* scale */
325 gij = wi * ai - di;
326 }
327 gij /= s;
328 gretl_matrix_set(dinfo->G, i, j, gij);
329 if (g != NULL) {
330 g[j] += gij;
331 }
332 }
333 }
334
335 return err;
336 }
337
338 #define matrix_plus(m,i,j,x) (m->val[(j)*m->rows+(i)]+=x)
339
340 /* Analytical Hessian: see Kalbfleisch and Prentice, 2002, pp. 69-70.
341 We're actually constructing the negative inverse of the Hessian here.
342 */
343
duration_hessian(double * theta,gretl_matrix * H,void * data)344 static int duration_hessian (double *theta,
345 gretl_matrix *H,
346 void *data)
347 {
348 duration_info *dinfo = (duration_info *) data;
349 const double *logt = dinfo->logt->val;
350 const double *Xb = dinfo->Xb->val;
351 double s, s2, Ai, wi, ewi, hwi;
352 int np = dinfo->npar;
353 double xij, xik, hjk;
354 int i, j, k, di;
355 int err = 0;
356
357 gretl_matrix_zero(H);
358
359 if (dinfo->dist == DUR_EXPON) {
360 s2 = s = 1;
361 } else {
362 s = theta[np - 1];
363 s2 = s * s;
364 }
365
366 for (i=0; i<dinfo->n; i++) {
367 di = uncensored(dinfo, i);
368 wi = (logt[i] - Xb[i]) / s;
369 ewi = exp(wi);
370
371 if (dinfo->dist == DUR_LOGLOG) {
372 Ai = (1 + di) * ewi / ((1 + ewi) * (1 + ewi));
373 } else if (dinfo->dist == DUR_LOGNORM) {
374 if (di) {
375 Ai = 1;
376 } else {
377 hwi = normal_h(wi);
378 Ai = hwi * (hwi - wi);
379 }
380 } else {
381 /* Weibull */
382 Ai = ewi;
383 }
384
385 for (j=0; j<np; j++) {
386 if (j < dinfo->k) {
387 /* covariate coeffs cross-block */
388 xij = gretl_matrix_get(dinfo->X, i, j);
389 for (k=0; k<=j; k++) {
390 xik = gretl_matrix_get(dinfo->X, i, k);
391 hjk = xij * xik * Ai / s2;
392 matrix_plus(H, j, k, hjk);
393 }
394 if (dinfo->dist != DUR_EXPON) {
395 /* coeff j and scale */
396 hjk = xij * wi * Ai / s2;
397 hjk += gretl_matrix_get(dinfo->G, i, j) / s;
398 matrix_plus(H, np - 1, j, hjk);
399 }
400 } else {
401 /* scale */
402 hjk = (wi * wi * Ai + di) / s2;
403 hjk += (2/s) * gretl_matrix_get(dinfo->G, i, j) / s;
404 matrix_plus(H, j, j, hjk);
405 }
406 }
407 }
408
409 /* fill out upper triangle */
410 gretl_matrix_mirror(H, 'L');
411
412 return err;
413 }
414
duration_hessian_inverse(double * theta,void * data,int * err)415 static gretl_matrix *duration_hessian_inverse (double *theta,
416 void *data,
417 int *err)
418 {
419 duration_info *dinfo = data;
420 gretl_matrix *H;
421
422 H = gretl_matrix_alloc(dinfo->npar, dinfo->npar);
423
424 if (H == NULL) {
425 *err = E_ALLOC;
426 } else {
427 *err = duration_hessian(theta, H, data);
428 }
429
430 #if 0
431 /* debugging check on numerical_hessian() */
432 gretl_matrix *nH = gretl_zero_matrix_new(dinfo->npar, dinfo->npar);
433
434 numerical_hessian(theta, nH, duration_loglik, data, 1, 0.0);
435 gretl_matrix_subtract_from(nH, H);
436 gretl_matrix_print(nH, "H: numerical - analytical");
437 gretl_matrix_free(nH);
438 #endif
439
440 if (!*err) {
441 *err = gretl_invert_symmetric_matrix(H);
442 }
443
444 return H;
445 }
446
447 /* calculate the OPG matrix at the starting point and use
448 its inverse (if any) as initial curvature matrix for BFGS
449 */
450
duration_init_H(duration_info * dinfo)451 static gretl_matrix *duration_init_H (duration_info *dinfo)
452 {
453 gretl_matrix *H = NULL;
454 int err;
455
456 dinfo->flags = DUR_UPDATE_XB;
457 err = duration_score(dinfo->theta, NULL, dinfo->npar,
458 NULL, dinfo);
459 dinfo->flags = 0;
460 if (!err) {
461 H = gretl_matrix_GG_inverse(dinfo->G, &err);
462 }
463
464 return H;
465 }
466
467 /* This is the last thing we do, after transcribing the
468 MLE results to pmod, so we don't have to worry about
469 saving and then restoring all the original values
470 attached to dinfo, which we will shortly destroy.
471 */
472
473 static void
duration_overall_LR_test(MODEL * pmod,duration_info * dinfo,int use_bfgs)474 duration_overall_LR_test (MODEL *pmod, duration_info *dinfo,
475 int use_bfgs)
476 {
477 double llu = dinfo->ll;
478 int err = 0;
479
480 dinfo->k = 1;
481 dinfo->npar = 1 + (dinfo->dist != DUR_EXPON);
482
483 gretl_matrix_reuse(dinfo->X, -1, dinfo->k);
484 gretl_matrix_reuse(dinfo->G, -1, dinfo->npar);
485 gretl_matrix_reuse(dinfo->beta, dinfo->k, 1);
486
487 dinfo->flags |= DUR_CONST_ONLY;
488 err = duration_estimates_init(dinfo);
489
490 /* now estimate constant-only model */
491
492 if (!err && use_bfgs) {
493 int maxit, fncount = 0, grcount = 0;
494 double toler;
495
496 BFGS_defaults(&maxit, &toler, DURATION);
497 err = BFGS_max(dinfo->theta, dinfo->npar, maxit, toler,
498 &fncount, &grcount, duration_loglik, C_LOGLIK,
499 duration_score, dinfo, NULL, OPT_NONE, NULL);
500 } else if (!err) {
501 double crittol = 1.0e-7;
502 double gradtol = 1.0e-7;
503 int iters = 0, maxit = 100;
504
505 err = newton_raphson_max(dinfo->theta, dinfo->npar, maxit,
506 crittol, gradtol, &iters,
507 C_LOGLIK, duration_loglik,
508 duration_score, duration_hessian,
509 dinfo, OPT_NONE, NULL);
510 }
511
512 if (!err && llu > dinfo->ll) {
513 pmod->chisq = 2 * (llu - dinfo->ll);
514 }
515 }
516
517 /*
518 Produce something to fill $yhat and $uhat for a duration model.
519
520 For $yhat we produce by default the conditional expectation E[t |
521 X, theta]. This carries the complication that the expected value
522 may be undefined for the log-logistic distribution. Alternatively,
523 if @opt includes OPT_M, we produce the conditional median, which is
524 always defined.
525
526 For $uhat we write Cox-Snell generalized residuals, namely the
527 integrated hazard function, which equals -log S(t).
528 */
529
duration_set_predictions(MODEL * pmod,duration_info * dinfo,const DATASET * dset,gretlopt opt)530 static void duration_set_predictions (MODEL *pmod, duration_info *dinfo,
531 const DATASET *dset, gretlopt opt)
532 {
533 const double *y = dset->Z[pmod->list[1]];
534 const double *logt = dinfo->logt->val;
535 int medians = (opt & OPT_M);
536 double St, G = 1.0;
537 double s = 1.0, p = 1.0;
538 double s22 = 0.0;
539 double pi_alpha = NADBL;
540 double wi, Xbi, expXbi;
541 int i, t;
542
543 if (dinfo->dist != DUR_EXPON) {
544 /* scale factor */
545 s = dinfo->theta[dinfo->npar-1];
546 p = 1 / s;
547 }
548
549 /* observation-invariant auxiliary quantities */
550
551 if (dinfo->dist == DUR_WEIBULL) {
552 /* agrees with Stata; R's "survreg" has this wrong? */
553 if (medians) {
554 G = pow(log(2.0), s);
555 } else {
556 G = gammafun(1 + s);
557 }
558 } else if (dinfo->dist == DUR_EXPON) {
559 if (medians) {
560 G = log(2.0);
561 } else {
562 G = gammafun(2.0);
563 }
564 } else if (dinfo->dist == DUR_LOGNORM) {
565 s22 = s * s / 2;
566 } else if (dinfo->dist == DUR_LOGLOG) {
567 if (!medians && s < 1) {
568 pi_alpha = M_PI * s / sin(M_PI * s);
569 }
570 }
571
572 i = 0;
573 for (t=pmod->t1; t<=pmod->t2; t++) {
574 if (na(pmod->yhat[t])) {
575 continue;
576 }
577
578 Xbi = dinfo->Xb->val[i];
579 wi = (logt[i] - Xbi) / s;
580 expXbi = exp(Xbi);
581
582 if (dinfo->dist == DUR_WEIBULL || dinfo->dist == DUR_EXPON) {
583 pmod->yhat[t] = expXbi * G;
584 St = exp(-exp(wi));
585 } else if (dinfo->dist == DUR_LOGNORM) {
586 if (medians) {
587 pmod->yhat[t] = expXbi;
588 } else {
589 pmod->yhat[t] = exp(Xbi + s22);
590 }
591 St = normal_cdf(-wi);
592 } else {
593 /* log-logistic */
594 if (medians) {
595 pmod->yhat[t] = expXbi;
596 } else if (s < 1) {
597 pmod->yhat[t] = expXbi * pi_alpha;
598 } else {
599 /* the expectation is undefined */
600 pmod->yhat[t] = NADBL;
601 }
602 St = 1.0 / (1 + pow(y[t] / expXbi, p));
603 }
604
605 /* generalized (Cox-Snell) residual */
606 pmod->uhat[t] = -log(St);
607 #if 0
608 if (!uncensored(dinfo, i)) {
609 pmod->uhat[t] += 1.0;
610 }
611 #endif
612
613 i++;
614 }
615
616 if (medians) {
617 pmod->opt |= OPT_M;
618 }
619 }
620
duration_model_add_vcv(MODEL * pmod,duration_info * dinfo,const DATASET * dset,gretlopt opt)621 static int duration_model_add_vcv (MODEL *pmod, duration_info *dinfo,
622 const DATASET *dset, gretlopt opt)
623 {
624 gretl_matrix *H = NULL;
625 int err = 0;
626
627 if (opt & OPT_G) {
628 err = gretl_model_add_OPG_vcv(pmod, dinfo->G, NULL);
629 } else {
630 H = duration_hessian_inverse(dinfo->theta, dinfo, &err);
631 if (!err) {
632 if (opt & OPT_R) {
633 err = gretl_model_add_QML_vcv(pmod, DURATION,
634 H, dinfo->G,
635 dset, opt, NULL);
636 } else {
637 err = gretl_model_add_hessian_vcv(pmod, H);
638 }
639 }
640 }
641
642 gretl_matrix_free(H);
643
644 return err;
645 }
646
647 static int
transcribe_duration_results(MODEL * pmod,duration_info * dinfo,const DATASET * dset,int fncount,int grcount,int use_bfgs,int censvar,gretlopt opt)648 transcribe_duration_results (MODEL *pmod, duration_info *dinfo,
649 const DATASET *dset,
650 int fncount, int grcount, int use_bfgs,
651 int censvar, gretlopt opt)
652 {
653 int np = dinfo->npar;
654 int j, v, err = 0;
655
656 pmod->ci = DURATION;
657
658 if (dinfo->dist == DUR_EXPON) {
659 pmod->opt |= OPT_E;
660 } else if (dinfo->dist == DUR_LOGLOG) {
661 pmod->opt |= OPT_L;
662 } else if (dinfo->dist == DUR_LOGNORM) {
663 pmod->opt |= OPT_Z;
664 } else {
665 pmod->opt |= OPT_B; /* Weibull */
666 }
667
668 if (censvar > 0) {
669 gretl_model_set_int(pmod, "cens_var", censvar);
670 }
671
672 if (use_bfgs) {
673 gretl_model_set_int(pmod, "fncount", fncount);
674 gretl_model_set_int(pmod, "grcount", grcount);
675 } else {
676 gretl_model_set_int(pmod, "iters", fncount);
677 }
678
679 if (!err) {
680 err = gretl_model_allocate_param_names(pmod, np);
681 if (!err) {
682 for (j=0; j<dinfo->k; j++) {
683 v = pmod->list[j+2];
684 gretl_model_set_param_name(pmod, j, dset->varname[v]);
685 }
686 if (dinfo->dist != DUR_EXPON) {
687 gretl_model_set_param_name(pmod, np-1, "sigma");
688 }
689 }
690 }
691
692 if (dinfo->dist == DUR_EXPON) {
693 pmod->sigma = 1.0;
694 } else {
695 pmod->sigma = dinfo->theta[np-1];
696 }
697
698 err = gretl_model_write_coeffs(pmod, dinfo->theta, np);
699
700 if (!err) {
701 err = duration_model_add_vcv(pmod, dinfo, dset, opt);
702 }
703
704 if (!err) {
705 duration_set_predictions(pmod, dinfo, dset, opt);
706 pmod->lnL = dinfo->ll;
707 mle_criteria(pmod, 0);
708 /* mask invalid statistics */
709 pmod->fstt = pmod->chisq = NADBL;
710 pmod->rsq = pmod->adjrsq = NADBL;
711 pmod->ess = NADBL;
712 /* but add overall LR test if possible */
713 if (np > 1 + (dinfo->dist != DUR_EXPON)) {
714 duration_overall_LR_test(pmod, dinfo, use_bfgs);
715 }
716 }
717
718 return err;
719 }
720
duration_estimate(MODEL * pmod,int censvar,const DATASET * dset,gretlopt opt,PRN * prn)721 int duration_estimate (MODEL *pmod, int censvar, const DATASET *dset,
722 gretlopt opt, PRN *prn)
723 {
724 gretlopt maxopt = (opt & OPT_V) | OPT_U;
725 duration_info dinfo;
726 int maxit = 200;
727 int fncount = 0;
728 int grcount = 0;
729 int use_bfgs = 0;
730 int err = 0;
731
732 if (opt & OPT_C) {
733 /* cluster implies robust */
734 opt |= OPT_R;
735 }
736
737 err = duration_init(&dinfo, pmod, censvar, dset, opt, prn);
738 if (err) {
739 goto bailout;
740 }
741
742 #if DDEBUG
743 if (!err && censvar > 0) {
744 fprintf(stderr, "duration: using var %d (%s) for censoring info\n",
745 censvar, dset->varname[censvar]);
746 }
747 #endif
748
749 if (libset_get_int(GRETL_OPTIM) == OPTIM_BFGS) {
750 use_bfgs = 1;
751 }
752
753 if (use_bfgs) {
754 /* initialize BFGS curvature */
755 gretl_matrix *H = duration_init_H(&dinfo);
756 double toler;
757
758 BFGS_defaults(&maxit, &toler, DURATION);
759 err = BFGS_max(dinfo.theta, dinfo.npar, maxit, toler,
760 &fncount, &grcount, duration_loglik, C_LOGLIK,
761 duration_score, &dinfo, H, maxopt,
762 dinfo.prn);
763 gretl_matrix_free(H);
764 } else {
765 double crittol = 1.0e-7;
766 double gradtol = 1.0e-7;
767
768 err = newton_raphson_max(dinfo.theta, dinfo.npar, maxit,
769 crittol, gradtol, &fncount,
770 C_LOGLIK, duration_loglik,
771 duration_score, duration_hessian,
772 &dinfo, maxopt, dinfo.prn);
773 }
774
775 if (!err) {
776 err = transcribe_duration_results(pmod, &dinfo, dset,
777 fncount, grcount, use_bfgs,
778 censvar, opt);
779 }
780
781 bailout:
782
783 duration_free(&dinfo);
784
785 if (err && !pmod->errcode) {
786 pmod->errcode = err;
787 }
788
789 return pmod->errcode;
790 }
791