1
2 #include "VARbase.h"
3 #include "VARio.h"
4 #include "dw_sbvar_command_line.h"
5 #include "sbvar_estimate.h"
6 #include "dw_switch.h"
7 #include "dw_switchio.h"
8 #include "dw_matrix_rand.h"
9 #include "dw_error.h"
10 #include "dw_ascii.h"
11 #include "dw_std.h"
12 #include "string.h"
13
SimulateData_VAR(TStateModel * model)14 TMatrix SimulateData_VAR(TStateModel *model)
15 {
16 T_VAR_Parameters *p=(T_VAR_Parameters*)(model->theta);
17 TVector initial, prob, *shocks;
18 TMatrix forecast;
19 int m, i, j;
20
21 // Allocate memory
22 prob=CreateVector(p->nstates);
23 initial=EquateVector((TVector)NULL,p->X[p->nobs]);
24 forecast=CreateMatrix(p->nobs,p->nvars);
25 shocks=dw_CreateArray_vector(p->nobs);
26 for (i=p->nobs-1; i >= 0; i--) shocks[i]=CreateVector(p->nvars);
27
28 // Get filtered probability at time T
29 for (j=p->nstates-1; j >= 0; j--)
30 ElementV(prob,j)=ProbabilityStateConditionalCurrent(j,p->nobs,model);
31
32 // Draw time T regime
33 p->S[0]=m=DrawDiscrete(prob);
34
35 // Draw regimes from time T+1 through T+h inclusive
36 for (j=0; j < p->nobs; j++)
37 {
38 ColumnVector(prob,model->sv->Q,m);
39 p->S[j]=m=DrawDiscrete(prob);
40 }
41
42 // Draw shocks
43 for (j=p->nobs-1; j >= 0; j--) dw_NormalVector(shocks[j]);
44
45 // Compute forecast
46 forecast_base(forecast,p->nobs,initial,shocks,p->S+1,model);
47
48 // Clean up and return
49 FreeVector(initial);
50 FreeVector(prob);
51 dw_FreeArray(shocks);
52
53 return forecast;
54 }
55
ActualData_VAR(TStateModel * model)56 TMatrix ActualData_VAR(TStateModel *model)
57 {
58 T_VAR_Parameters *p=(T_VAR_Parameters*)(model->theta);
59 TMatrix Y;
60 int i, j;
61
62 // Allocate memory
63 Y=CreateMatrix(p->nobs,p->nvars);
64
65 for (i=p->nobs-1; i >= 0; i--)
66 for (j=p->nvars-1; j >= 0; j--)
67 ElementM(Y,i,j)=ElementV(p->Y[i+1],j);
68
69 return Y;
70 }
71
main(int nargs,char ** args)72 int main(int nargs, char **args)
73 {
74 TStateModel *model, *model_constant;
75 TVARCommandLine *cmd=(TVARCommandLine*)NULL, *cmd_constant=(TVARCommandLine*)NULL;
76 TEstimateInfo* info;
77 T_VAR_Parameters *p, *theta;
78 TMatrix X, Y, A0, Aplus;
79 int **translation_table;
80 FILE *f_out;
81 char *filename, *header;
82 int i, j, s;
83
84 // Setup model
85 if (!(cmd=CreateTStateModel_VARCommandLine(nargs,args,(TVARCommandLine*)NULL)))
86 dw_Error(MEM_ERR);
87 else
88 if (!cmd->model)
89 dw_UserError("Unable to setup SBVAR from command line parameters");
90 else
91 {
92 p=(T_VAR_Parameters*)(cmd->model->theta);
93
94 if (p->npre - p->nvars*p->nlags > 1)
95 printf("Only constant is allowed as an exogenous variable\n");
96 else
97 {
98 // Set new out tag
99 if (!cmd->out_tag)
100 if (!cmd->in_tag)
101 cmd->out_tag=dw_DuplicateString("constant_notag");
102 else
103 {
104 cmd->out_tag=(char*)dw_malloc((strlen(cmd->in_tag)+10)*sizeof(char));
105 strcat(strcpy(cmd->out_tag,"constant_"),cmd->in_tag);
106 }
107
108 filename=CreateFilenameFromTag("%sinit_%s.dat",cmd->out_tag,cmd->in_directory);
109 if (f_out=fopen(filename,"wt"))
110 {
111 // Simulate Data
112 X=CreateMatrix(p->nobs,p->npre);
113 Y=SimulateData_VAR(cmd->model);
114 for (j=p->npre-1; j >= 0; j--) ElementM(X,0,j)=ElementV(p->X[p->nobs],j);
115
116 // The below lines will use the actual data
117 // Y=ActualData_VAR(cmd->model);
118 // for (j=p->npre-1; j >= 0; j--) ElementM(X,0,j)=ElementV(p->X[1],j);
119
120 for (i=1; i < p->nobs; i++)
121 {
122 for (j=p->nvars*(p->nlags-1)-1; j >= 0; j--) ElementM(X,i,j+p->nvars)=ElementM(X,i-1,j);
123 for (j=p->nvars-1; j >= 0; j--) ElementM(X,i,j)=ElementM(Y,i-1,j);
124 for (j=p->nvars*p->nlags; j < p->npre; j++) ElementM(X,i,j)=ElementM(X,i-1,j);
125 }
126
127 // Create new TStateModel
128 theta=CreateTheta_VAR(p->Specification,p->nvars,p->nlags,p->npre - p->nvars*p->nlags,p->nstates,p->nobs,p->coef_states,p->var_states,p->U,p->V,p->W,Y,X);
129 if (p->Specification & SPEC_SIMS_ZHA)
130 SetPriors_VAR_SimsZha(theta,p->A0_prior,p->Aplus_prior,p->zeta_a_prior,p->zeta_b_prior,p->lambda_prior);
131 else
132 SetPriors_VAR(theta,p->A0_prior,p->Aplus_prior,p->zeta_a_prior,p->zeta_b_prior);
133 model=CreateStateModel(p->nobs,DuplicateMarkovStateVariable(cmd->model->sv),CreateRoutines_VAR(),theta);
134
135 // Write specification
136 Write_VAR_Specification(f_out,(char*)NULL,model);
137
138 // Write true parameters
139 header="Truth: ";
140 fprintf(f_out,"//== %sFree parameters ==//\n",header);
141 WriteFreeParameters(f_out,(char*)NULL,cmd->model);
142 fprintf(f_out,"\n");
143 WriteBaseTransitionMatrices(f_out,(char*)NULL,header,cmd->model);
144 Write_VAR_Parameters(f_out,(char*)NULL,header,cmd->model);
145
146 // Create constant parameter model and estimate
147 dw_InitializeArray_int(translation_table=dw_CreateRectangularArray_int(p->nvars,1),0);
148 theta=CreateTheta_VAR(p->Specification,p->nvars,p->nlags,p->npre - p->nvars*p->nlags,1,p->nobs,translation_table,translation_table,p->U,p->V,p->W,Y,X);
149 dw_FreeArray(translation_table);
150 if (p->Specification & SPEC_SIMS_ZHA)
151 SetPriors_VAR_SimsZha(theta,p->A0_prior,p->Aplus_prior,p->zeta_a_prior,p->zeta_b_prior,p->lambda_prior);
152 else
153 SetPriors_VAR(theta,p->A0_prior,p->Aplus_prior,p->zeta_a_prior,p->zeta_b_prior);
154 model_constant=CreateStateModel(p->nobs,CreateMarkovStateVariable_ConstantState(0),CreateRoutines_VAR(),theta);
155 InitializeParameters_VAR((T_VAR_Parameters*)model_constant->theta);
156 ThetaChanged(model_constant);
157 cmd_constant=Base_VARCommandLine(nargs,args);
158 cmd_constant->model=model_constant;
159 if (cmd_constant->out_tag) dw_free(cmd_constant->out_tag);
160 cmd_constant->out_tag=dw_DuplicateString("constant_simulated_data");
161 info=SetupEstimateFromVARCommandLine(nargs,args,cmd_constant);
162 FindMode_VAR_csminwel(info,(FILE*)NULL);
163
164 // Insert constant parameters
165 p=(T_VAR_Parameters*)(model->theta);
166 // A0
167 A0=MakeA0((TMatrix)NULL,0,(T_VAR_Parameters*)(model_constant->theta));
168 for (j=p->nvars-1; j >= 0; j--)
169 for (s=p->n_coef_states[j]-1; s >= 0; s--)
170 for (i=p->nvars-1; i >= 0; i--)
171 ElementV(p->A0[j][s],i)=ElementM(A0,i,j);
172 FreeMatrix(A0);
173 // Aplus
174 Aplus=MakeAplus((TMatrix)NULL,0,(T_VAR_Parameters*)(model_constant->theta));
175 for (j=p->nvars-1; j >= 0; j--)
176 for (s=p->n_coef_states[j]-1; s >= 0; s--)
177 for (i=p->npre-1; i >= 0; i--)
178 ElementV(p->Aplus[j][s],i)=ElementM(Aplus,i,j);
179 FreeMatrix(Aplus);
180 // Zeta
181 for (j=p->nvars-1; j >= 0; j--)
182 for (s=p->n_var_states[j]-1; s >= 0; s--)
183 p->Zeta[j][s]=1.0;
184 // b0, bplus, lambda, and psi
185 Update_b0_bplus_from_A0_Aplus(p);
186 if (p->Specification & SPEC_SIMS_ZHA) Update_lambda_psi_from_bplus(p);
187 // Flags
188 p->valid_parameters=1;
189 // Transition matrix
190 DefaultTransitionMatrixParameters(model);
191
192 // Write initial values
193 header="Initial: ";
194 WriteBaseTransitionMatrices(f_out,(char*)NULL,header,model);
195 Write_VAR_Parameters(f_out,(char*)NULL,header,model);
196
197 // Clean-up
198 FreeMatrix(X);
199 FreeMatrix(Y);
200 FreeStateModel(model);
201 FreeEstimateInfo(info);
202 }
203 dw_free(filename);
204 }
205 }
206
207 Free_VARCommandLine(cmd);
208 return 0;
209 }
210