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_matrix.h"
19 #include "dw_rand.h"
20 #include "dw_parse_cmd.h"
21 #include "dw_ascii.h"
22 #include "dw_error.h"
23 #include "VARbase.h"
24 #include "VARio.h"
25 #include "switch.h"
26 #include "switchio.h"
27 #include "dw_sbvar_command_line.h"
28 #include "dw_std.h"
29
30 #include <stdlib.h>
31 #include <string.h>
32 #include <math.h>
33 #include <time.h>
34
35 /*
36 Returns the largest of the absolute values of the eigenvalues of the companion
37 matrix of B
38 */
LargestRoot(TMatrix B)39 static PRECISION LargestRoot(TMatrix B)
40 {
41 TMatrix X;
42 TVector re_eig, im_eig;
43 int i, k=RowM(B), n=ColM(B);
44 PRECISION max, tmp;
45 InitializeMatrix(X=CreateMatrix(k,k),0.0);
46 for (i=n; i < k; i++) ElementM(X,i-n,i)=1.0;
47 InsertSubMatrix(X,B,0,0,0,0,k,n);
48 re_eig=CreateVector(k);
49 im_eig=CreateVector(k);
50 Eigenvalues(re_eig,im_eig,X);
51 max=ElementV(re_eig,0)*ElementV(re_eig,0)+ElementV(im_eig,0)*ElementV(im_eig,0);
52 for (i=k-1; i > 0; i--)
53 if ((tmp=ElementV(re_eig,i)*ElementV(re_eig,i)+ElementV(im_eig,i)*ElementV(im_eig,i)) > max) max=tmp;;
54 FreeVector(im_eig);
55 FreeVector(re_eig);
56 FreeMatrix(X);
57 return sqrt(max);
58 }
59
60 /*
61 Computes the matrices A0(s)*Xi(s) and B(s) = Ahat(s) * inverse(A0(s))
62 where Ahat(s) are the first nvars*nlags rows of Aplus(s), for 0 <= s < h.
63 The integer h is the number of regimes. Both A0_Xi and B must be matrix
64 arrays of length h and each element of A0_Xi and B must be of the proper
65 size.
66 */
GetIRMatrices(int h,TMatrix * A0_Xi,TMatrix * B,T_VAR_Parameters * p)67 void GetIRMatrices(int h, TMatrix *A0_Xi, TMatrix *B, T_VAR_Parameters *p)
68 {
69 TMatrix Aplus;
70 int i, j, k;
71 PRECISION tmp;
72 for (k=h-1; k >= 0; k--)
73 {
74 MakeA0(A0_Xi[k],k,p);
75 Aplus=MakeAplus((TMatrix)NULL,k,p);
76 SubMatrix(B[k],Aplus,0,0,p->nvars*p->nlags,p->nvars);
77 FreeMatrix(Aplus);
78 ProductInverseMM(B[k],B[k],A0_Xi[k]);
79 for (j=p->nvars-1; j >= 0; j--)
80 for (tmp=sqrt(p->Zeta[j][p->var_states[j][k]]), i=p->nvars-1; i >= 0; i--)
81 ElementM(A0_Xi[k],i,j)*=tmp;
82 }
83 }
84
85
86 /*
87 Additional command line parameters
88
89 -ndraws <integer>
90 Number of draws to save (default value = 1000)
91 -burnin <integer>
92 Number of burn-in draws (default = 0.1 * thin * ndraws)
93 -thin <integer>
94 Thinning factor. Total number of draws made is thin*ndraws + burnin (default = 1)
95 -no_deg <integer : retry> <integer : cutoff>
96 Rejects degenerate draws of regimes. Will redraw regimes up to retry times. In
97 order not to be degenerate, there must be at least cutoff observations for each
98 state. (default = 1000 1)
99 -append
100 Append simulation and info to existing files if they exist and use last parameter
101 draw as initial value. (default overwrites existing files)
102 -mh <integer>
103 Tuning period for Metropolis-Hasting draws (default value = 30000)
104 -flat
105 Produce flat output file and header file (default value = off)
106 -nofree
107 Do not produce free parameters file (default value = off)
108 -nd1
109 Normalize diagonal of A0 to one (flat output only)
110 */
111 //#define DEBUG_1
dw_sbvar_simulate_command_line(int nargs,char ** args)112 void dw_sbvar_simulate_command_line(int nargs, char **args)
113 {
114 TStateModel *model;
115 T_VAR_Parameters *p;
116 FILE *f_flat, *f_free, *f_info;
117 char *filename;
118 time_t begin_time, end_time;
119 int count, check, period=1000, i, tuning, burnin, ndraws, thin, append, writeflat, writefree, nd1;
120 TVARCommandLine *cmd=(TVARCommandLine*)NULL;
121 TVector parameters;
122
123 // Get tuning peroid, burn-in period, number of iterations, and thinning factor
124 ndraws=dw_ParseInteger_String(nargs,args,"ndraws",1000);
125 thin=dw_ParseInteger_String(nargs,args,"thin",1);
126 burnin=dw_ParseInteger_String(nargs,args,"burnin",(ndraws*thin)/10);
127 tuning=dw_ParseInteger_String(nargs,args,"mh",30000);
128
129 // Append file
130 append=(dw_FindArgument_String(nargs,args,"append") >= 0) ? 1 : 0;
131
132 // Get output options
133 writeflat=(dw_FindArgument_String(nargs,args,"flat") >= 0) ? 1 : 0;
134 writefree=(dw_FindArgument_String(nargs,args,"nofree") >= 0) ? 0 : 1;
135 nd1=(dw_FindArgument_String(nargs,args,"nd1") >= 0) ? 1 : 0;
136
137 // Setup model and initial parameters
138 if (!(cmd=CreateTStateModel_VARCommandLine(nargs,args,(TVARCommandLine*)NULL)))
139 dw_Error(MEM_ERR);
140 else
141 if (!cmd->model)
142 dw_UserError("Unable to setup SBVAR from command line parameters");
143 else
144 if (!cmd->parameters_filename_actual)
145 dw_UserError("Unable to initialize parameters from command line parameters");
146 else
147 {
148 model=cmd->model;
149 p=(T_VAR_Parameters*)(model->theta);
150
151 // Degeneracy behavior
152 if ((i=dw_FindArgument_String(nargs,args,"no_deg")) == -1)
153 model->reject_degenerate_draws=0;
154 else
155 {
156 if ((i+1 >= nargs) || ((model->reject_degenerate_draws=atoi(args[i+1])) <= 0)) model->reject_degenerate_draws=100;
157 if ((i+2 >= nargs) || ((model->degenerate_draws_cutoff=atoi(args[i+2])) <= 0)) model->degenerate_draws_cutoff=1;
158 }
159
160
161 // Allocate parameters workspace
162 parameters=CreateVector(NumberFreeParametersTheta(model)+NumberFreeParametersQ(model));
163
164 // initialize parameters from previous simulation run if append
165 if (append)
166 if (writefree)
167 {
168 filename=CreateFilenameFromTag("%ssimulation_%s.out",cmd->out_tag,cmd->out_directory);
169 f_free=fopen(filename,"rt");
170 dw_free(filename);
171 if (f_free)
172 {
173 while (GetPosteriorDraw(f_free,model,F_FREE));
174 fclose(f_free);
175 }
176 }
177 else
178 if (writeflat)
179 {
180 filename=CreateFilenameFromTag("%sdraws_%s.out",cmd->out_tag,cmd->out_directory);
181 f_flat=fopen(filename,"rt");
182 dw_free(filename);
183 if (f_flat)
184 {
185 while (GetPosteriorDraw(f_flat,model,F_FLAT));
186 fclose(f_flat);
187 }
188 }
189
190 // Save initial parameters
191 ConvertThetaToFreeParameters(model,pElementV(parameters));
192 ConvertQToFreeParameters(model,pElementV(parameters)+NumberFreeParametersTheta(model));
193
194 // Output tag
195 if (!cmd->out_tag)
196 cmd->out_tag=dw_DuplicateString(cmd->in_tag ? cmd->in_tag : "");
197
198 // Open header and flat file and print headers
199 if (writeflat)
200 {
201 filename=CreateFilenameFromTag("%sdraws_header_%s.out",cmd->out_tag,cmd->out_directory);
202 if (f_flat=fopen(filename,"wt"))
203 {
204 dw_free(filename);
205 WriteBaseTransitionMatricesFlat_Headers(f_flat,model);
206 Write_VAR_ParametersFlat_Headers(f_flat,model);
207 fprintf(f_flat,"\n");
208 fclose(f_flat);
209 }
210
211 filename=CreateFilenameFromTag("%sdraws_%s.out",cmd->out_tag,cmd->out_directory);
212 f_flat=fopen(filename,append ? "at" : "wt");
213 dw_free(filename);
214 }
215 else
216 f_flat=(FILE*)NULL;
217
218 // Open free file
219 if (writefree)
220 {
221 filename=CreateFilenameFromTag("%ssimulation_%s.out",cmd->out_tag,cmd->out_directory);
222 f_free=fopen(filename,append ? "at" : "wt");
223 dw_free(filename);
224 }
225 else
226 f_free=(FILE*)NULL;
227
228 // Begin time
229 begin_time=(int)time((time_t*)NULL);
230
231 // Open info file and write
232 filename=CreateFilenameFromTag("%ssimulation_info_%s.out",cmd->out_tag,cmd->out_directory);
233 f_info=fopen(filename,append ? "at" : "wt");
234 dw_free(filename);
235 if (f_info)
236 {
237 fprintf(f_info,"\nSimulation Run Info\n");
238 fprintf(f_info,"Number saved draws: %d\n",ndraws);
239 fprintf(f_info,"Total draws: %d\n",ndraws*thin);
240 fprintf(f_info,"Thinning factor: %d\n",thin);
241 fprintf(f_info,"Burnin period: %d\n",burnin);
242 fprintf(f_info,"Adaptive tuning period: %d\n",tuning);
243 fprintf(f_info,"File format: %s\n",f_free ? (f_flat ? "free and flat formats" : "free format") : (f_flat ? "flat format" : "file not written"));
244 fprintf(f_info,"Seed: %d\n",dw_ParseInteger_String(nargs,args,"seed",0));
245 fprintf(f_info,"Begin time stamp: %s",ctime(&begin_time));
246 }
247
248 // Calibration of jumping parameters
249 printf("Calibrating jumping parameters - %d draws\n",tuning);
250 AdaptiveMetropolisScale(model,tuning,1000,1,(FILE*)NULL); // tuning iterations - 1000 iterations before updating - verbose
251 end_time=(int)time((time_t*)NULL);
252 printf("Elapsed Time: %d seconds\n",(int)(end_time - begin_time));
253
254 // Reset initial parameters
255 ConvertFreeParametersToTheta(model,pElementV(parameters));
256 ConvertFreeParametersToQ(model,pElementV(parameters)+NumberFreeParametersTheta(model));
257
258 // Reset degeneracy counts
259 model->n_degenerate_draws=model->total_degenerate_draws=0;
260
261 // Burn-in period
262 printf("Burn-in period - %d draws\n",burnin);
263 for (check=period, count=1; count <= burnin; count++)
264 {
265 DrawAll(model);
266
267 if (count == check)
268 {
269 check+=period;
270 printf("%d iterations completed out of %d\n",count,burnin);
271 end_time=(int)time((time_t*)NULL);
272 printf("Elapsed Time: %d seconds\n",(int)(end_time - begin_time));
273 }
274 }
275 end_time=(int)time((time_t*)NULL);
276 printf("Elapsed Time: %d seconds\n",(int)(end_time - begin_time));
277 ResetMetropolisInformation(p);
278
279 #ifdef DEBUG_1
280 TMatrix *A0_Xi, *B;
281 FILE *f_debug;
282
283 f_debug=fopen("debug.out","wt");
284 //f_debug=stdout;
285
286 A0_Xi=dw_CreateArray_matrix(model->sv->nstates);
287 B=dw_CreateArray_matrix(model->sv->nstates);
288 for (i=model->sv->nstates-1; i >= 0; i--)
289 {
290 A0_Xi[i]=CreateMatrix(p->nvars,p->nvars);
291 B[i]=CreateMatrix(p->nvars*p->nlags,p->nvars);
292 }
293 fprintf(f_debug,"%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model));
294 /* GetIRMatrices(model->sv->nstates,A0_Xi,B,p); */
295 /* for (i=0; i < model->sv->nstates; i++) */
296 /* fprintf(f_debug," %lf",LargestRoot(B[i])); */
297 /* fprintf(f_debug,"\n"); */
298 DrawStates(model);
299 #endif
300
301 // Simulation
302 printf("Simulating - %d draws\n",ndraws);
303 for (check=period, count=1; count <= ndraws; count++)
304 {
305 #ifdef DEBUG_1
306 if (f_free)
307 {
308 //printf("Previous: %lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model));
309
310 //DrawStates(model);
311 //printf("States: %lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model));
312
313 /* DrawTransitionMatrixParameters(model); */
314 /* printf("%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
315
316 /* DrawTheta(model); */
317 /* printf("%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
318
319 /* DrawZeta_DotProducts(model); */
320 /* ((T_VAR_Parameters*)(model->theta))->valid_parameters=1; */
321 /* ThetaChanged(model); */
322 /* printf("Zeta: %lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
323
324 /* DrawA0_Metropolis(model); */
325 /* ((T_VAR_Parameters*)(model->theta))->valid_parameters=1; */
326 /* ThetaChanged(model); */
327 /* printf("A0: %lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
328
329 /* DrawAplus(model); */
330 /* ((T_VAR_Parameters*)(model->theta))->valid_parameters=1; */
331 /* ThetaChanged(model); */
332 /* printf("%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
333
334 if (p->Specification & SPEC_SIMS_ZHA)
335 {
336 Draw_psi(model);
337 Update_bplus_from_lambda_psi(p);
338 Update_Aplus_from_bplus_A0(p);
339 p->valid_parameters=1;
340 ThetaChanged(model);
341 fprintf(f_debug,"%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model));
342
343 /* Draw_lambda(model); */
344 /* Update_bplus_from_lambda_psi(p); */
345 /* Update_Aplus_from_bplus_A0(p); */
346 /* p->valid_parameters=1; */
347 /* ThetaChanged(model); */
348 /* printf("lambda: %lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
349 }
350 else
351 { printf("Not Sim-Zha specification\n"); dw_exit(0); }
352
353 // Normalize
354 //Normalize_VAR((T_VAR_Parameters*)(model->theta));
355
356 //getchar();
357
358 /* fprintf(f_free,"%lf %lf %lf\n",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model),LogPrior(model)); */
359 /* GetIRMatrices(model->sv->nstates,A0_Xi,B,p); */
360 /* for (i=0; i < model->sv->nstates; i++) */
361 /* fprintf(f_free," %lf",LargestRoot(B[i])); */
362 /* fprintf(f_free,"\n"); */
363 }
364 #else
365 for (i=thin; i > 0; i--) DrawAll(model);
366
367 if (f_flat)
368 {
369 WriteBaseTransitionMatricesFlat(f_flat,(char*)NULL,model,"%lf ");
370 if (nd1)
371 Write_VAR_ParametersFlat_A0_Diagonal_One(f_flat,model,"%lf ");
372 else
373 Write_VAR_ParametersFlat(f_flat,model,"%lf ");
374 fprintf(f_flat,"\n");
375 }
376
377 if (f_free)
378 {
379 fprintf(f_free,"%le %le ",LogPosterior_StatesIntegratedOut(model),LogLikelihood_StatesIntegratedOut(model));
380 ConvertThetaToFreeParameters(model,pElementV(parameters));
381 ConvertQToFreeParameters(model,pElementV(parameters)+NumberFreeParametersTheta(model));
382 dw_PrintVector(f_free,parameters,"%le ");
383 }
384 #endif
385
386 if (count == check)
387 {
388 check+=period;
389 printf("%d * %d iterations completed out of %d * %d\n",count,thin,ndraws,thin);
390 end_time=(int)time((time_t*)NULL);
391 printf("Elapsed Time: %d seconds\n",(int)(end_time - begin_time));
392 }
393 }
394 end_time=(int)time((time_t*)NULL);
395 printf("Total Elapsed Time: %d seconds\n",(int)(end_time - begin_time));
396
397 if (f_info)
398 {
399 fprintf(f_info,"End time stamp: %s\n",ctime(&end_time));
400 fprintf(f_info,"Number degenerate draws of regimes: %d (%d)\n",model->n_degenerate_draws,model->total_degenerate_draws);
401 fprintf(f_info,"degenerate draws %s\n",model->reject_degenerate_draws ? "rejected" : "kept");
402 /* fprintf(f_info,"Degenerate draws by regime: "); */
403 /* for (i=0; i < model->sv->nstates; i++) fprintf(f_info,"%d ",model->states_count[i]); */
404 /* fprintf(f_info,"\n"); */
405 }
406
407 // clean up
408 if (f_flat) fclose(f_flat);
409 if (f_free) fclose(f_free);
410 if (f_info) fclose(f_info);
411 FreeVector(parameters);
412 }
413
414 // clean up
415 Free_VARCommandLine(cmd);
416 }
417