1 /*
2 * Copyright (C) 1996-2011 Daniel Waggoner and Tao Zha
3 *
4 * This free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * It is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * If you did not received a copy of the GNU General Public License
15 * with this software, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "dw_state_space_impulse_response.h"
19 #include "dw_MSStateSpace.h"
20 #include "dw_metropolis_simulation.h"
21 #include "dw_matrix_array.h"
22 #include "dw_histogram.h"
23 #include "dw_parse_cmd.h"
24 #include "dw_math.h"
25 #include "dw_std.h"
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <string.h>
30
31 /*
32 Assumes:
33 IR : (h+1) x (nz*nepsilon) matrix
34 S : integer array of length h+1
35 h : non-negative integer
36 model : valid pointer to TStateModel/T_MSStateSpace structure
37
38 Returns:
39 One upon success and zero otherwise.
40
41 Results:
42 Fills IR with the impulse responses of the state variables with respect to
43 the fundamental shocks for horizons 0 through h inclusive.
44
45 Notes:
46 The current value of the parameters implicit in model are used to form the
47 impulse responses. The element in position (k,i+j*nz) is the response of
48 the ith state variable to the jth fundamental shock at horizon k. Horizon 0
49 is the contemporaneous response.
50 */
51 /* int dw_state_space_impulse_response_states(TMatrix IR, int *S, int h, TStateModel *model) */
52 /* { */
53 /* TVector *z, e; */
54 /* T_MSStateSpace *statespace=(T_MSStateSpace*)(model->theta); */
55 /* int i, j, k; */
56
57 /* // initialize structure */
58 /* if ((statespace->t0 < 0) && !Initialize_MSStateSpace(model)) return 0; */
59
60 /* // check that there are fundamental shocks */
61 /* if (statespace->nepsilon == 0) return 0; */
62
63 /* // allocate memory */
64 /* z=dw_CreateArray_vector(statespace->nepsilon); */
65
66 /* // contemporaneous response to fundamental shocks */
67 /* InitializeVector(e=CreateVector(statespace->nepsilon),0.0); */
68 /* for (j=statespace->nepsilon-1; j >= 0; j--) */
69 /* { */
70 /* ElementV(e,j)=1.0; */
71 /* z[j]=ProductMV((TVector)NULL,statespace->Phiz[S[0]],e); */
72 /* for (i=statespace->nz-1; i >= 0; i--) */
73 /* ElementM(IR,0,i+j*statespace->nz)=ElementV(z[j],i); */
74 /* ElementV(e,j)=0.0; */
75 /* } */
76 /* FreeVector(e); */
77
78 /* // horizon k response to fundamental shocks */
79 /* for (k=1; k <= h; k++) */
80 /* for (j=statespace->nepsilon-1; j >= 0; j--) */
81 /* { */
82 /* ProductMV(z[j],statespace->F[S[k]],z[j]); */
83 /* for (i=statespace->nz-1; i >= 0; i--) */
84 /* ElementM(IR,k,i+j*statespace->nz)=ElementV(z[j],i); */
85 /* } */
86
87 /* // free memory and return success */
88 /* dw_FreeArray(z); */
89 /* return 1; */
90 /* } */
91
92 /*
93 Assumes:
94 IR : (h+1) x (ny*(nepsilon+nu)) matrix
95 S : integer array of length h+1
96 h : non-negative integer
97 model : valid pointer to TStateModel/T_MSStateSpace structure
98
99 Returns:
100 One upon success and zero otherwise.
101
102 Results:
103 Fills IR with the impulse responses of the observable variables with respect
104 to the fundamental and measurement error shocks for horizons 0 through h
105 inclusive.
106
107 Notes:
108 The current value of the parameters implicit in model are used to form the
109 impulse responses. The element in position (k,i+j*ny) is the response of
110 the ith observable variable to the jth shock at horizon k. Shocks 0 through
111 nepsilon-1 are fundamental and shocks nepsilon through nepsilon+nu-1 are
112 measurement error. Horizon 0 is the contemporaneous response.
113 */
114 /* int dw_state_space_impulse_response_observables(TMatrix IR, int *S, int h, TStateModel *model) */
115 /* { */
116 /* TVector *z, y, e; */
117 /* T_MSStateSpace *statespace=(T_MSStateSpace*)(model->theta); */
118 /* int i, j, k; */
119
120 /* // initialize structure */
121 /* if ((statespace->t0 < 0) && !Initialize_MSStateSpace(model)) return 0; */
122
123 /* // check that there are fundamental or measurement shocks */
124 /* if (statespace->nepsilon + statespace->nu == 0) return 0; */
125
126 /* // allocate memory and initialize */
127 /* y=CreateVector(statespace->ny); */
128 /* InitializeMatrix(IR,0.0); */
129
130 /* // contemporaneous response to fundamental shocks */
131 /* if (statespace->nepsilon > 0) */
132 /* { */
133 /* z=dw_CreateArray_vector(statespace->nepsilon); */
134 /* InitializeVector(e=CreateVector(statespace->nepsilon),0.0); */
135 /* for (j=statespace->nepsilon-1; j >= 0; j--) */
136 /* { */
137 /* ElementV(e,j)=1.0; */
138 /* z[j]=ProductMV((TVector)NULL,statespace->Phiz[S[0]],e); */
139 /* ProductMV(y,statespace->H[S[0]],z[j]); */
140 /* for (i=statespace->ny-1; i >= 0; i--) */
141 /* ElementM(IR,0,i+j*statespace->ny)=ElementV(y,i); */
142 /* ElementV(e,j)=0.0; */
143 /* } */
144 /* FreeVector(e); */
145 /* } */
146
147 /* // contemporaneous response to measument error shocks */
148 /* if (statespace->nu > 0) */
149 /* { */
150 /* InitializeVector(e=CreateVector(statespace->nu),0.0); */
151 /* for (j=statespace->nu-1; j >= 0; j--) */
152 /* { */
153 /* ElementV(e,j)=1.0; */
154 /* ProductMV(y,statespace->Phiy[S[0]],e); */
155 /* for (i=statespace->ny-1; i >= 0; i--) */
156 /* ElementM(IR,0,i+(statespace->nepsilon+j)*statespace->ny)=ElementV(y,i); */
157 /* ElementV(e,j)=0.0; */
158 /* } */
159 /* FreeVector(e); */
160 /* } */
161
162 /* // horizon k response to fundamental shocks */
163 /* if (statespace->nepsilon > 0) */
164 /* { */
165 /* for (k=1; k <= h; k++) */
166 /* for (j=statespace->nepsilon-1; j >= 0; j--) */
167 /* { */
168 /* ProductMV(z[j],statespace->F[S[k]],z[j]); */
169 /* ProductMV(y,statespace->H[S[k]],z[j]); */
170 /* for (i=statespace->ny-1; i >= 0; i--) */
171 /* ElementM(IR,k,i+j*statespace->ny)=ElementV(y,i); */
172 /* } */
173 /* dw_FreeArray(z); */
174 /* } */
175
176 /* // free memory and return success */
177 /* FreeVector(y); */
178 /* return 1; */
179 /* } */
180
181 /*
182 Assumes
183 f_out : valid FILE pointer
184 s : constant base regime for the impluse response.
185 percentiles : vector of numbers between 0 and 1.
186 posterior_draws : each row contains log posterior, log likelihood, and posterior draw or null
187 thin : thinning factor
188 h : non-negative integer
189 model : point to valid TStateModel/T_MSStateSpace structure
190
191 Results:
192 Computes and prints to the file f_out the requested percentiles of the
193 impulse responses of the state variables to the fundamental shocks given
194 one is always in regime s.
195
196 Returns:
197 One upon success and zero otherwise.
198
199 Notes:
200 For each percentile, the element in position (k,i+j*nz) is the response of
201 the ith state variable to the jth shock at horizon k. Shocks 0 through
202 nepsilon-1 are fundamental and shocks nepsilon through nepsilon+nu-1 are
203 measurement error. Horizon 0 is the contemporaneous response.
204 */
dw_state_space_impulse_response_states_percentile_regime(FILE * f_out,int s,TVector percentiles,TMatrix posterior_draws,int thin,int h,TStateModel * model)205 int dw_state_space_impulse_response_states_percentile_regime(FILE *f_out, int s, TVector percentiles, TMatrix posterior_draws,
206 int thin, int h, TStateModel *model)
207 {
208 T_MSStateSpace *statespace;
209 int rtrn=0, ntheta, nq, *S, i;
210 TVector x;
211 TMatrix IR;
212 TMatrixHistogram *histogram;
213
214 // quick check of passed parameters
215 if (!f_out || !percentiles || (thin <= 0) || (h < 0) || !model) return 0;
216
217 statespace=(T_MSStateSpace*)(model->theta);
218 ntheta=NumberFreeParametersTheta(model);
219 nq=NumberFreeParametersQ(model);
220
221 // s must be in the proper range
222 if ((s < 0) || (statespace->nbasestates <= s)) return 0;
223
224 // allocate memory
225 S=(int*)dw_malloc((h+1)*sizeof(int));
226 IR=CreateMatrix(h+1,statespace->nz*statespace->nepsilon);
227 histogram=CreateMatrixHistogram(RowM(IR),ColM(IR),40,HISTOGRAM_VARIABLE);
228
229 // set S to s
230 for (i=0; i <= h; i++) S[i]=s;
231
232 if (posterior_draws)
233 {
234 x=CreateVector(ColM(posterior_draws));
235 for (i=RowM(posterior_draws)-1; i >= 0; i-=thin)
236 {
237 // Get free parameters and push them into model
238 RowVector(x,posterior_draws,i);
239 if (ntheta > 0) ConvertFreeParametersToTheta(model,pElementV(x)+2);
240 if (nq > 0) ConvertFreeParametersToQ(model,pElementV(x)+2+ntheta);
241
242 // Get transition matrix -- this should be modified if transition matrix is time varying
243 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
244 goto ERROR_EXIT;
245
246 // Compute impulse response
247 if (!dw_state_space_impulse_response(IR,S,model,IR_STATES))
248 goto ERROR_EXIT;
249
250 // Accumulate impulse response
251 AddMatrixObservation(IR,histogram);
252 }
253 FreeVector(x);
254 }
255 else
256 {
257 // Get transition matrix -- this should be modified if transition matrix is time varying
258 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
259 goto ERROR_EXIT;
260
261 // Compute impulse response
262 if (!dw_state_space_impulse_response(IR,S,model,IR_STATES))
263 goto ERROR_EXIT;
264
265 // Accumulate impulse response
266 AddMatrixObservation(IR,histogram);
267 }
268
269 for (i=0; i < DimV(percentiles); i++)
270 {
271 MatrixPercentile(IR,ElementV(percentiles,i),histogram);
272 dw_PrintMatrix(f_out,IR,"%lg ");
273 fprintf(f_out,"\n");
274 }
275
276 rtrn=1;
277
278 ERROR_EXIT:
279 FreeMatrixHistogram(histogram);
280 FreeMatrix(IR);
281 dw_free(S);
282
283 return rtrn;
284 }
285
286 /*
287 Assumes
288 f_out : valid FILE pointer
289 s : constant regime for the impluse response.
290 percentiles : vector of numbers between 0 and 1 inclusive
291 posterior_draws : each row contains log posterior, log likelihood, and posterior draw or null
292 thin : thinning factor
293 h : non-negative integer
294 model : point to valid TStateModel/T_MSStateSpace structure
295
296 Results:
297 Computes and prints to the file f_out the requested percentiles of the
298 impulse responses of the observables variables to the fundamental and
299 measurement error shocks.
300
301 Returns:
302 One upon success and zero otherwise.
303
304 Notes:
305 For each percentile, the element in position (k,i+j*ny) is the response of
306 the ith observable variable to the jth shock at horizon k. Shocks 0 through
307 nepsilon-1 are fundamental and shocks nepsilon through nepsilon+nu-1 are
308 measurement error. Horizon 0 is the contemporaneous response.
309 */
dw_state_space_impulse_response_observables_percentile_regime(FILE * f_out,int s,TVector percentiles,TMatrix posterior_draws,int thin,int h,TStateModel * model)310 int dw_state_space_impulse_response_observables_percentile_regime(FILE *f_out, int s, TVector percentiles, TMatrix posterior_draws,
311 int thin, int h, TStateModel *model)
312 {
313 T_MSStateSpace *statespace;
314 int rtrn=0, ntheta, nq, *S, i;
315 TVector x;
316 TMatrix IR;
317 TMatrixHistogram *histogram;
318
319 // quick check of passed parameters
320 if (!f_out || !percentiles || (thin <= 0) || (h < 0) || !model) return 0;
321
322 statespace=(T_MSStateSpace*)(model->theta);
323 ntheta=NumberFreeParametersTheta(model);
324 nq=NumberFreeParametersQ(model);
325
326 // s must be in the proper range
327 if ((s < 0) || (statespace->nbasestates <= s)) return 0;
328
329 // allocate memory
330 S=(int*)dw_malloc((h+1)*sizeof(int));
331 IR=CreateMatrix(h+1,statespace->ny*(statespace->nepsilon+statespace->nu));
332 histogram=CreateMatrixHistogram(RowM(IR),ColM(IR),40,HISTOGRAM_VARIABLE);
333
334 // set S to s
335 for (i=0; i <= h; i++) S[i]=s;
336
337 if (posterior_draws)
338 {
339 x=CreateVector(ColM(posterior_draws));
340 for (i=RowM(posterior_draws)-1; i >= 0; i-=thin)
341 {
342 // Get free parameters and push them into model
343 RowVector(x,posterior_draws,i);
344 if (ntheta > 0) ConvertFreeParametersToTheta(model,pElementV(x)+2);
345 if (nq > 0) ConvertFreeParametersToQ(model,pElementV(x)+2+ntheta);
346
347 // Get transition matrix -- this should be modified if transition matrix is time varying
348 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
349 goto ERROR_EXIT;
350
351 // Compute impulse response
352 if (!dw_state_space_impulse_response(IR,S,model,IR_OBSERVABLES))
353 goto ERROR_EXIT;
354
355 // Accumulate impulse response
356 AddMatrixObservation(IR,histogram);
357 }
358 FreeVector(x);
359 }
360 else
361 {
362 // Get transition matrix -- this should be modified if transition matrix is time varying
363 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
364 goto ERROR_EXIT;
365
366 // Compute impulse response
367 if (!dw_state_space_impulse_response(IR,S,model,IR_OBSERVABLES))
368 goto ERROR_EXIT;
369
370 // Accumulate impulse response
371 AddMatrixObservation(IR,histogram);
372 }
373
374 for (i=0; i < DimV(percentiles); i++)
375 {
376 MatrixPercentile(IR,ElementV(percentiles,i),histogram);
377 dw_PrintMatrix(f_out,IR,"%lg ");
378 fprintf(f_out,"\n");
379 }
380
381 rtrn=1;
382
383 ERROR_EXIT:
384 FreeMatrixHistogram(histogram);
385 FreeMatrix(IR);
386 dw_free(S);
387
388
389 return rtrn;
390 }
391
392 /*
393 Assumes
394 f_out : valid FILE pointer
395 percentiles : vector of numbers between 0 and 1 inclusive
396 posterior_draws : each row contains log posterior, log likelihood, and posterior draw
397 thin : thinning factor
398 init_prob : vector or null pointer
399 h : non-negative integer
400 model : point to valid TStateModel/T_MSStateSpace structure
401
402 Results:
403 Computes and prints to the file f_out the requested percentiles of the
404 impulse responses of the state variables to the fundamental shocks.
405
406 Returns:
407 One upon success and zero otherwise.
408
409 Notes:
410 For each percentile, the element in position (k,i+j*nz) is the response of
411 the ith state variable to the jth shock at horizon k. Shocks 0 through
412 nepsilon-1 are fundamental and shocks nepsilon through nepsilon+nu-1 are
413 measurement error. Horizon 0 is the contemporaneous response.
414 */
dw_state_space_impulse_response_states_percentile_ergodic(FILE * f_out,TVector percentiles,TMatrix posterior_draws,int thin,TVector init_prob,int h,TStateModel * model)415 int dw_state_space_impulse_response_states_percentile_ergodic(FILE *f_out, TVector percentiles, TMatrix posterior_draws,
416 int thin, TVector init_prob, int h, TStateModel *model)
417 {
418 T_MSStateSpace *statespace;
419 int rtrn=0, ntheta, nq, *S, use_ergodic, i, j, m;
420 TVector x, prob;
421 TMatrix IR;
422 TMatrixHistogram *histogram;
423
424 // quick check of passed parameters
425 if (!f_out || !percentiles || !posterior_draws || (thin <= 0) || (h < 0) || !model) return 0;
426
427 statespace=(T_MSStateSpace*)(model->theta);
428 ntheta=NumberFreeParametersTheta(model);
429 nq=NumberFreeParametersQ(model);
430
431 // allocate memory
432 x=CreateVector(ColM(posterior_draws));
433 prob=CreateVector(model->sv->nstates);
434 S=(int*)dw_malloc((h+1)*sizeof(int));
435 IR=CreateMatrix(h+1,statespace->nz*statespace->nepsilon);
436 histogram=CreateMatrixHistogram(RowM(IR),ColM(IR),40,HISTOGRAM_VARIABLE);
437
438 // Use ergodic initial probability
439 if (!init_prob)
440 {
441 use_ergodic=1;
442 init_prob=CreateVector(model->sv->nstates);
443 }
444 else
445 use_ergodic=0;
446
447 for (i=RowM(posterior_draws)-1; i >= 0; i-=thin)
448 {
449 // Get free parameters and push them into model
450 RowVector(x,posterior_draws,i);
451 if (ntheta > 0) ConvertFreeParametersToTheta(model,pElementV(x)+2);
452 if (nq > 0) ConvertFreeParametersToQ(model,pElementV(x)+2+ntheta);
453
454 // Get transition matrix -- this should be modified if transition matrix is time varying
455 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
456 goto ERROR_EXIT;
457
458 // Get ergodic distribution if necessary
459 if (use_ergodic) Ergodic(init_prob,model->sv->Q);
460
461 // Draw regimes
462 m=DrawDiscrete(init_prob);
463 S[0]=model->sv->lag_index[m][0];
464 for (j=1; j <= h; j++)
465 {
466 ColumnVector(prob,model->sv->Q,m);
467 m=DrawDiscrete(prob);
468 S[j]=model->sv->lag_index[m][0];
469 }
470
471 // Compute impulse response
472 if (!dw_state_space_impulse_response(IR,S,model,IR_STATES))
473 goto ERROR_EXIT;
474
475 // Accumulate impulse response
476 AddMatrixObservation(IR,histogram);
477 }
478
479 for (i=0; i < DimV(percentiles); i++)
480 {
481 MatrixPercentile(IR,ElementV(percentiles,i),histogram);
482 dw_PrintMatrix(f_out,IR,"%lg ");
483 fprintf(f_out,"\n");
484 }
485
486 rtrn=1;
487
488 ERROR_EXIT:
489 if (use_ergodic) dw_free(init_prob);
490 FreeMatrixHistogram(histogram);
491 FreeMatrix(IR);
492 dw_free(S);
493 FreeVector(prob);
494 FreeVector(x);
495
496 return rtrn;
497 }
498
499 /*
500 Assumes
501 f_out : valid FILE pointer
502 percentiles : vector of numbers between 0 and 1 inclusive
503 posterior_draws : each row contains log posterior, log likelihood, and posterior draw
504 thin : thinning factor
505 init_prob : vector or null pointer
506 h : non-negative integer
507 model : point to valid TStateModel/T_MSStateSpace structure
508
509 Results:
510 Computes and prints to the file f_out the requested percentiles of the
511 impulse responses of the observables variables to the fundamental and
512 measurement error shocks.
513
514 Returns:
515 One upon success and zero otherwise.
516
517 Notes:
518 For each percentile, the element in position (k,i+j*ny) is the response of
519 the ith observable variable to the jth shock at horizon k. Shocks 0 through
520 nepsilon-1 are fundamental and shocks nepsilon through nepsilon+nu-1 are
521 measurement error. Horizon 0 is the contemporaneous response.
522 */
dw_state_space_impulse_response_observables_percentile_ergodic(FILE * f_out,TVector percentiles,TMatrix posterior_draws,int thin,TVector init_prob,int h,TStateModel * model)523 int dw_state_space_impulse_response_observables_percentile_ergodic(FILE *f_out, TVector percentiles, TMatrix posterior_draws,
524 int thin, TVector init_prob, int h, TStateModel *model)
525 {
526 T_MSStateSpace *statespace;
527 int rtrn=0, ntheta, nq, *S, use_ergodic, i, j, m;
528 TVector x, prob;
529 TMatrix IR;
530 TMatrixHistogram *histogram;
531
532 // quick check of passed parameters
533 if (!f_out || !percentiles || !posterior_draws || (thin <= 0) || (h < 0) || !model) return 0;
534
535 statespace=(T_MSStateSpace*)(model->theta);
536 ntheta=NumberFreeParametersTheta(model);
537 nq=NumberFreeParametersQ(model);
538
539 // allocate memory
540 x=CreateVector(ColM(posterior_draws));
541 prob=CreateVector(model->sv->nstates);
542 S=(int*)dw_malloc((h+1)*sizeof(int));
543 IR=CreateMatrix(h+1,statespace->ny*(statespace->nepsilon+statespace->nu));
544 histogram=CreateMatrixHistogram(RowM(IR),ColM(IR),40,HISTOGRAM_VARIABLE);
545
546 // Use ergodic initial probability
547 if (!init_prob)
548 {
549 use_ergodic=1;
550 init_prob=CreateVector(model->sv->nstates);
551 }
552 else
553 use_ergodic=0;
554
555 for (i=RowM(posterior_draws)-1; i >= 0; i-=thin)
556 {
557 // Get free parameters and push them into model
558 RowVector(x,posterior_draws,i);
559 if (ntheta > 0) ConvertFreeParametersToTheta(model,pElementV(x)+2);
560 if (nq > 0) ConvertFreeParametersToQ(model,pElementV(x)+2+ntheta);
561
562 // Get transition matrix -- this should be modified if transition matrix is time varying
563 if ((model->sv->t1 < model->nobs) && !sv_ComputeTransitionMatrix(model->nobs,model->sv,model))
564 goto ERROR_EXIT;
565
566 // Get ergodic distribution if necessary
567 if (use_ergodic) Ergodic(init_prob,model->sv->Q);
568
569 // Draw regimes
570 m=DrawDiscrete(init_prob);
571 S[0]=model->sv->lag_index[m][0];
572 for (j=1; j <= h; j++)
573 {
574 ColumnVector(prob,model->sv->Q,m);
575 m=DrawDiscrete(prob);
576 S[j]=model->sv->lag_index[m][0];
577 }
578
579 // Compute impulse response
580 if (!dw_state_space_impulse_response(IR,S,model,IR_OBSERVABLES))
581 goto ERROR_EXIT;
582
583 // Accumulate impulse response
584 AddMatrixObservation(IR,histogram);
585 }
586
587 for (i=0; i < DimV(percentiles); i++)
588 {
589 MatrixPercentile(IR,ElementV(percentiles,i),histogram);
590 dw_PrintMatrix(f_out,IR,"%lg ");
591 fprintf(f_out,"\n");
592 }
593
594 rtrn=1;
595
596 ERROR_EXIT:
597 if (use_ergodic) dw_free(init_prob);
598 FreeMatrixHistogram(histogram);
599 FreeMatrix(IR);
600 dw_free(S);
601 FreeVector(prob);
602 FreeVector(x);
603
604 return rtrn;
605 }
606
607