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 "gretl_matrix.h"
23
24 #include "gretl_f2c.h"
25 #include "clapack_double.h"
26
27 #include <gtk/gtk.h>
28
29 struct flag_info {
30 GtkWidget *dialog;
31 GtkWidget *levcheck;
32 GtkWidget *infcheck;
33 GtkWidget *dffcheck;
34 unsigned char *flag;
35 };
36
destroy_save_dialog(GtkWidget * w,struct flag_info * finfo)37 static gboolean destroy_save_dialog (GtkWidget *w, struct flag_info *finfo)
38 {
39 free(finfo);
40 gtk_main_quit();
41 return FALSE;
42 }
43
update_save_flag(GtkWidget * w,struct flag_info * finfo)44 static gboolean update_save_flag (GtkWidget *w, struct flag_info *finfo)
45 {
46 int checked = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(w));
47
48 if (w == finfo->levcheck) {
49 if (checked) {
50 *finfo->flag |= SAVE_LEVERAGE;
51 } else {
52 *finfo->flag &= ~SAVE_LEVERAGE;
53 }
54 } else if (w == finfo->infcheck) {
55 if (checked) {
56 *finfo->flag |= SAVE_INFLUENCE;
57 } else {
58 *finfo->flag &= ~SAVE_INFLUENCE;
59 }
60 } else if (w == finfo->dffcheck) {
61 if (checked) {
62 *finfo->flag |= SAVE_DFFITS;
63 } else {
64 *finfo->flag &= ~SAVE_DFFITS;
65 }
66 }
67
68 return FALSE;
69 }
70
cancel_set_flag(GtkWidget * w,struct flag_info * finfo)71 static gboolean cancel_set_flag (GtkWidget *w, struct flag_info *finfo)
72 {
73 *finfo->flag = 0;
74 gtk_widget_destroy(finfo->dialog);
75 return FALSE;
76 }
77
save_dialog_finalize(GtkWidget * w,struct flag_info * finfo)78 static gboolean save_dialog_finalize (GtkWidget *w, struct flag_info *finfo)
79 {
80 gtk_widget_destroy(finfo->dialog);
81 return FALSE;
82 }
83
leverage_data_dialog(void)84 unsigned char leverage_data_dialog (void)
85 {
86 struct flag_info *finfo;
87 GtkWidget *dialog, *tmp, *button, *vbox, *hbox;
88 GtkWidget *internal_vbox;
89 unsigned char flag = SAVE_LEVERAGE | SAVE_INFLUENCE | SAVE_DFFITS;
90
91 finfo = malloc(sizeof *finfo);
92 if (finfo == NULL) return 0;
93
94 dialog = gtk_dialog_new();
95
96 finfo->dialog = dialog;
97 finfo->flag = &flag;
98
99 vbox = gtk_dialog_get_content_area(GTK_DIALOG(dialog));
100 hbox = gtk_dialog_get_action_area(GTK_DIALOG(dialog));
101
102 gtk_window_set_title(GTK_WINDOW(dialog), _("gretl: save data"));
103 gtk_window_set_resizable(GTK_WINDOW(dialog), FALSE);
104 gtk_container_set_border_width(GTK_CONTAINER(vbox), 10);
105 gtk_container_set_border_width(GTK_CONTAINER(hbox), 5);
106 gtk_box_set_spacing(GTK_BOX(vbox), 5);
107
108 gtk_window_set_position(GTK_WINDOW(dialog), GTK_WIN_POS_MOUSE);
109
110 g_signal_connect(G_OBJECT(dialog), "destroy",
111 G_CALLBACK(destroy_save_dialog), finfo);
112
113 internal_vbox = gtk_vbox_new(FALSE, 5);
114
115 hbox = gtk_hbox_new(FALSE, 5);
116 tmp = gtk_label_new(_("Variables to save:"));
117 gtk_box_pack_start(GTK_BOX(hbox), tmp, TRUE, TRUE, 5);
118 gtk_box_pack_start(GTK_BOX(internal_vbox), hbox, TRUE, TRUE, 5);
119
120 /* Leverage */
121 button = gtk_check_button_new_with_label(_("leverage"));
122 gtk_box_pack_start(GTK_BOX(internal_vbox), button, TRUE, TRUE, 0);
123 gtk_toggle_button_set_active (GTK_TOGGLE_BUTTON (button), TRUE);
124 g_signal_connect(G_OBJECT(button), "clicked",
125 G_CALLBACK(update_save_flag), finfo);
126 finfo->levcheck = button;
127
128 /* Influence */
129 button = gtk_check_button_new_with_label(_("influence"));
130 gtk_box_pack_start(GTK_BOX(internal_vbox), button, TRUE, TRUE, 0);
131 gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON (button), TRUE);
132 g_signal_connect(G_OBJECT(button), "clicked",
133 G_CALLBACK(update_save_flag), finfo);
134 finfo->infcheck = button;
135
136 /* DFFITS */
137 button = gtk_check_button_new_with_label(_("DFFITS"));
138 gtk_box_pack_start(GTK_BOX(internal_vbox), button, TRUE, TRUE, 0);
139 gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON (button), TRUE);
140 g_signal_connect(G_OBJECT(button), "clicked",
141 G_CALLBACK(update_save_flag), finfo);
142 finfo->dffcheck = button;
143
144 hbox = gtk_hbox_new(FALSE, 5);
145 gtk_box_pack_start(GTK_BOX(hbox), internal_vbox, TRUE, TRUE, 5);
146
147 gtk_widget_show(internal_vbox);
148 gtk_box_pack_start(GTK_BOX(vbox), hbox, TRUE, TRUE, 5);
149
150 hbox = gtk_dialog_get_action_area(GTK_DIALOG(dialog));
151 gtk_button_box_set_layout(GTK_BUTTON_BOX(hbox), GTK_BUTTONBOX_END);
152 gtk_box_set_spacing(GTK_BOX(hbox), 10);
153
154 /* Cancel button */
155 tmp = gtk_button_new_from_stock(GTK_STOCK_CANCEL);
156 gtk_container_add(GTK_CONTAINER(hbox), tmp);
157 g_signal_connect(G_OBJECT(tmp), "clicked",
158 G_CALLBACK(cancel_set_flag), finfo);
159
160 /* "OK" button */
161 tmp = gtk_button_new_from_stock(GTK_STOCK_OK);
162 gtk_container_add(GTK_CONTAINER(hbox), tmp);
163 g_signal_connect(G_OBJECT(tmp), "clicked",
164 G_CALLBACK(save_dialog_finalize), finfo);
165 gtk_widget_set_can_default(tmp, TRUE);
166 gtk_widget_grab_default(tmp);
167
168 gtk_widget_show_all(dialog);
169
170 gtk_main();
171
172 return flag;
173 }
174
175 static void
leverage_x_range(int t1,int t2,const double * x,FILE * fp)176 leverage_x_range (int t1, int t2, const double *x, FILE *fp)
177 {
178 double xrange, xmin;
179 double xmin0 = x[t1];
180 double xmax = x[t2];
181
182 xrange = xmax - xmin0;
183 xmin = xmin0 - xrange * .025;
184 xmax += xrange * .025;
185
186 if (xmin0 >= 0.0 && xmin < 0.0) {
187 xmin = 0.0;
188 }
189
190 fprintf(fp, "set xrange [%.7g:%.7g]\n", xmin, xmax);
191 }
192
leverage_plot(const MODEL * pmod,gretl_matrix * S,DATASET * dset)193 static int leverage_plot (const MODEL *pmod, gretl_matrix *S,
194 DATASET *dset)
195 {
196 FILE *fp;
197 const double *obs = NULL;
198 int t, err = 0;
199
200 fp = open_plot_input_file(PLOT_LEVERAGE, 0, &err);
201 if (err) {
202 return err;
203 }
204
205 if (dataset_is_time_series(dset)) {
206 obs = gretl_plotx(dset, OPT_NONE);
207 if (obs == NULL) {
208 if (fp != NULL) {
209 fclose(fp);
210 }
211 return 1;
212 }
213 }
214
215 gretl_push_c_numeric_locale();
216
217 fputs("set size 1.0,1.0\nset multiplot\nset size 1.0,0.48\n", fp);
218 fputs("set xzeroaxis\n", fp);
219 fputs("set nokey\n", fp);
220
221 if (obs != NULL) {
222 leverage_x_range(pmod->t1, pmod->t2, obs, fp);
223 } else {
224 fprintf(fp, "set xrange [%g:%g]\n",
225 pmod->t1 + 0.5, pmod->t2 + 1.5);
226 }
227
228 /* upper plot: leverage factor */
229 fputs("set origin 0.0,0.50\n", fp);
230 gnuplot_missval_string(fp);
231 fputs("set yrange [0:1]\n", fp);
232 fprintf(fp, "set title '%s'\n", _("leverage"));
233 fputs("plot \\\n'-' using 1:2 w impulses\n", fp);
234
235 for (t=pmod->t1; t<=pmod->t2; t++) {
236 double h = gretl_matrix_get(S, t - pmod->t1, 0);
237
238 if (na(h)) {
239 if (obs != NULL) {
240 fprintf(fp, "%g ?\n", obs[t]);
241 } else {
242 fprintf(fp, "%d ?\n", t + 1);
243 }
244 } else {
245 if (obs != NULL) {
246 fprintf(fp, "%g %g\n", obs[t], h);
247 } else {
248 fprintf(fp, "%d %g\n", t + 1, h);
249 }
250 }
251 }
252 fputs("e\n", fp);
253
254 /* lower plot: influence factor */
255 fputs("set origin 0.0,0.0\n", fp);
256 gnuplot_missval_string(fp);
257 fputs("set yrange [*:*]\n", fp);
258 fprintf(fp, "set title '%s'\n", _("influence"));
259 fputs("plot \\\n'-' using 1:2 w impulses\n", fp);
260
261 for (t=pmod->t1; t<=pmod->t2; t++) {
262 double f = gretl_matrix_get(S, t - pmod->t1, 1);
263
264 if (na(f)) {
265 if (obs != NULL) {
266 fprintf(fp, "%g ?\n", obs[t]);
267 } else {
268 fprintf(fp, "%d ?\n", t + 1);
269 }
270 } else {
271 if (obs != NULL) {
272 fprintf(fp, "%g %g\n", obs[t], f);
273 } else {
274 fprintf(fp, "%d %g\n", t + 1, f);
275 }
276 }
277 }
278 fputs("e\n", fp);
279
280 fputs("unset multiplot\n", fp);
281
282 gretl_pop_c_numeric_locale();
283
284 return finalize_plot_input_file(fp);
285 }
286
leverage_print(const MODEL * pmod,gretl_matrix * S,double Xvalcrit,DATASET * dset,PRN * prn)287 static void leverage_print (const MODEL *pmod,
288 gretl_matrix *S,
289 double Xvalcrit,
290 DATASET *dset,
291 PRN *prn)
292 {
293 double lp = 2.0 * pmod->ncoeff / pmod->nobs;
294 int obslen = max_obs_marker_length(dset);
295 int t, j, gotlp = 0;
296
297 if (obslen < 8) {
298 obslen = 8;
299 }
300
301 bufspace(obslen, prn);
302 pprintf(prn, "%*s", UTF_WIDTH(_("residual"), 16), _("residual"));
303 pprintf(prn, "%*s", UTF_WIDTH(_("leverage"), 16), _("leverage"));
304 pprintf(prn, "%*s", UTF_WIDTH(_("influence"), 16), _("influence"));
305 pprintf(prn, "%*s", UTF_WIDTH(_("DFFITS"), 14), _("DFFITS"));
306 pputc(prn, '\n');
307 bufspace(obslen, prn);
308 pputs(prn, " u 0<=h<=1 u*h/(1-h)\n\n");
309
310 for (t=pmod->t1, j=0; t<=pmod->t2; t++, j++) {
311 double h, st, d, f;
312 gchar *fstr = NULL;
313
314 if (na(pmod->uhat[t])) {
315 print_obs_marker(t, dset, obslen, prn);
316 pputc(prn, '\n');
317 continue;
318 }
319
320 h = gretl_matrix_get(S, j, 0);
321 if (h > lp) {
322 gotlp = 1;
323 }
324
325 f = gretl_matrix_get(S, j, 1);
326 if (na(f)) {
327 fstr = g_strdup(_("undefined"));
328 gretl_utf8_truncate(fstr, 15);
329 } else {
330 fstr = g_strdup_printf("%15.5g", f);
331 }
332
333 print_obs_marker(t, dset, obslen, prn);
334
335 st = gretl_matrix_get(S, j, 2);
336 d = st * sqrt(h / (1.0 - h));
337 pprintf(prn, "%14.5g %14.3f%s %s %14.3f\n", pmod->uhat[t], h,
338 (h > lp)? "*" : " ", fstr, d);
339 g_free(fstr);
340 }
341
342 if (gotlp) {
343 pprintf(prn, "\n%s\n", _("('*' indicates a leverage point)"));
344 } else {
345 pprintf(prn, "\n%s\n", _("No leverage points were found"));
346 }
347
348 pprintf(prn, "\n%s = %g\n\n", _("Cross-validation criterion"), Xvalcrit);
349 }
350
model_leverage(const MODEL * pmod,DATASET * dset,gretlopt opt,PRN * prn,int * err)351 gretl_matrix *model_leverage (const MODEL *pmod, DATASET *dset,
352 gretlopt opt, PRN *prn,
353 int *err)
354 {
355 integer info, lwork;
356 integer m, n, lda;
357 gretl_matrix *Q, *S = NULL;
358 double *tau, *work;
359 double Xvalcrit, df_adj = 0;
360 int xvc_only = 0;
361 int i, j, s, t, vi;
362 /* allow for missing obs in model range */
363 int modn = pmod->t2 - pmod->t1 + 1;
364
365 m = pmod->nobs; /* # of rows = # of observations */
366 lda = m; /* leading dimension of Q */
367 n = pmod->list[0] - 1; /* # of cols = # of variables */
368
369 Q = gretl_matrix_alloc(m, n);
370
371 /* dim of tau is min (m, n) */
372 tau = malloc(n * sizeof *tau);
373 work = malloc(sizeof *work);
374
375 if (Q == NULL || tau == NULL || work == NULL) {
376 *err = E_ALLOC;
377 goto qr_cleanup;
378 }
379
380 /* copy independent var values into Q, skipping missing obs */
381 j = 0;
382 for (i=2; i<=pmod->list[0]; i++) {
383 vi = pmod->list[i];
384 for (t=pmod->t1; t<=pmod->t2; t++) {
385 if (!na(pmod->uhat[t])) {
386 Q->val[j++] = dset->Z[vi][t];
387 }
388 }
389 }
390
391 /* do a workspace size query */
392 lwork = -1;
393 info = 0;
394 dgeqrf_(&m, &n, Q->val, &lda, tau, work, &lwork, &info);
395 if (info != 0) {
396 fprintf(stderr, "dgeqrf: info = %d\n", (int) info);
397 *err = 1;
398 goto qr_cleanup;
399 }
400
401 /* set up optimally sized work array */
402 lwork = (integer) work[0];
403 work = realloc(work, (size_t) lwork * sizeof *work);
404 if (work == NULL) {
405 *err = E_ALLOC;
406 goto qr_cleanup;
407 }
408
409 /* run actual QR factorization */
410 dgeqrf_(&m, &n, Q->val, &lda, tau, work, &lwork, &info);
411 if (info != 0) {
412 fprintf(stderr, "dgeqrf: info = %d\n", (int) info);
413 *err = 1;
414 goto qr_cleanup;
415 }
416
417 /* obtain the real "Q" matrix */
418 dorgqr_(&m, &n, &n, Q->val, &lda, tau, work, &lwork, &info);
419 if (info != 0) {
420 *err = 1;
421 goto qr_cleanup;
422 }
423
424 /* If we're neither printing anything nor saving leverage
425 at al as series, then we just need the cross-validation
426 criterion.
427 */
428 if ((opt & OPT_Q) && !(opt & OPT_S)) {
429 xvc_only = 1; /* the S matrix is not needed */
430 } else {
431 S = gretl_matrix_alloc(modn, 3);
432 if (S == NULL) {
433 *err = E_ALLOC;
434 goto qr_cleanup;
435 }
436 gretl_matrix_set_t1(S, pmod->t1);
437 df_adj = sqrt(pmod->dfd - 1.0);
438 }
439
440 /* initialize cross-validation criterion */
441 Xvalcrit = 0.0;
442
443 /* do the "h" calculations, etc. */
444 s = 0;
445 for (t=pmod->t1, i=0; t<=pmod->t2; t++, i++) {
446 double den, f = NADBL, d = NADBL;
447 double q, h, et = pmod->uhat[t];
448
449 if (na(et)) {
450 h = NADBL;
451 } else {
452 h = 0.0;
453 for (j=0; j<n; j++) {
454 q = gretl_matrix_get(Q, s, j);
455 h += q * q;
456 }
457 if (h < 1.0) {
458 f = et / (1 - h);
459 Xvalcrit += f * f;
460 f -= et;
461 if (!xvc_only) {
462 /* studentized residual */
463 den = sqrt((1 - h) * pmod->ess - et * et);
464 d = df_adj / den * et;
465 }
466 }
467 s++;
468 }
469 if (S != NULL) {
470 gretl_matrix_set(S, i, 0, h);
471 gretl_matrix_set(S, i, 1, f);
472 gretl_matrix_set(S, i, 2, d);
473 }
474 }
475
476 record_test_result(Xvalcrit, NADBL);
477
478 /* print the results, unless in quiet mode */
479 if (!(opt & OPT_Q)) {
480 leverage_print(pmod, S, Xvalcrit, dset, prn);
481 if (gnuplot_graph_wanted(PLOT_LEVERAGE, opt)) {
482 leverage_plot(pmod, S, dset);
483 }
484 }
485
486 qr_cleanup:
487
488 gretl_matrix_free(Q);
489 free(tau);
490 free(work);
491
492 return S;
493 }
494