diff --git a/matlab/swz/c-code/sbvar/switching/switchio.c b/matlab/swz/c-code/sbvar/switching/switchio.c index 3f58083eb..7f246f45d 100644 --- a/matlab/swz/c-code/sbvar/switching/switchio.c +++ b/matlab/swz/c-code/sbvar/switching/switchio.c @@ -562,7 +562,7 @@ int ReadBaseTransitionMatrices_SV(FILE *f_in, TMarkovStateVariable* sv, char *he ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded); /* Free Q */ - FreeMatrix(Q); + FreeMatrix(Q); /* // Update ansi-c*/ if (!Update_B_from_Q_SV(sv)) @@ -728,42 +728,42 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) if (sv->n_state_variables > 1) { for (i=0; i < sv->n_state_variables; i++) - if (!ReadBaseTransitionMatricesFlat_SV(f_in,sv->state_variable[i])) - return 0; + if (!ReadBaseTransitionMatricesFlat_SV(f_in,sv->state_variable[i])) + return 0; } else { // Read transition matrix Q=CreateMatrix(sv->nbasestates,sv->nbasestates); for (j=0; j < ColM(Q); j++) - for (i=0; i < RowM(Q); i++) - if (fscanf(f_in," %lf ",&ElementM(Q,i,j)) != 1) - { - FreeMatrix(Q); - return 0; - } + for (i=0; i < RowM(Q); i++) + if (fscanf(f_in," %lf ",&ElementM(Q,i,j)) != 1) + { + FreeMatrix(Q); + return 0; + } // Scale the columns of Q - loose requirement on sumation to one for (j=sv->nbasestates-1; j >= 0; j--) - { - for (sum=0.0, i=sv->nbasestates-1; i >= 0; i--) - if (ElementM(Q,i,j) < 0.0) - { - FreeMatrix(Q); - dw_UserError("Transition matrix can not have negative elements."); - return 0; - } - else - sum+=ElementM(Q,i,j); - if (fabs(sum-1.0) > 1.0e-4) - { - FreeMatrix(Q); - dw_UserError("Transition matrix columns must sum to one."); - return 0; - } - for (sum=1.0/sum, i=sv->nbasestates-1; i >= 0; i--) - ElementM(Q,i,j)*=sum; - } + { + for (sum=0.0, i=sv->nbasestates-1; i >= 0; i--) + if (ElementM(Q,i,j) < 0.0) + { + FreeMatrix(Q); + dw_UserError("Transition matrix can not have negative elements."); + return 0; + } + else + sum+=ElementM(Q,i,j); + if (fabs(sum-1.0) > 1.0e-4) + { + FreeMatrix(Q); + dw_UserError("Transition matrix columns must sum to one."); + return 0; + } + for (sum=1.0/sum, i=sv->nbasestates-1; i >= 0; i--) + ElementM(Q,i,j)*=sum; + } // Convert base transition matrix to full transition matrix. ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded); @@ -773,10 +773,10 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) // Update if (!Update_B_from_Q_SV(sv)) - { - dw_UserError("Transition matrices do not satisfy restrictions"); - return 0; - } + { + dw_UserError("Transition matrices do not satisfy restrictions"); + return 0; + } return 1; } diff --git a/matlab/swz/c-code/sbvar/var/VARbase.c b/matlab/swz/c-code/sbvar/var/VARbase.c index bbed86e84..186ef20c3 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.c +++ b/matlab/swz/c-code/sbvar/var/VARbase.c @@ -2433,13 +2433,13 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s { ProductVM(y,x,Aplus[S[t]]); if (shocks) - { - for (i=p->nvars-1; i >= 0; i--) - ElementV(y,i)+=ElementV(shocks[t],i)/sqrt(p->Zeta[i][p->var_states[i][S[t]]]); - } + { + for (i=p->nvars-1; i >= 0; i--) + ElementV(y,i)+=ElementV(shocks[t],i)/sqrt(p->Zeta[i][p->var_states[i][S[t]]]); + } ProductInverseVM(y,y,A0[S[t]]); for (i=p->nvars-1; i >= 0; i--) - ElementM(forecast,t,i)=ElementV(y,i); + ElementM(forecast,t,i)=ElementV(y,i); memmove(pElementV(x)+p->nvars,pElementV(x),(p->nlags-1)*p->nvars*sizeof(PRECISION)); memcpy(pElementV(x),pElementV(y),p->nvars*sizeof(PRECISION)); @@ -2487,23 +2487,23 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s /* Pxi1[k]=CreateVector(statespace->nstates); */ /* SPxi[k]=CreateVector(statespace->nstates); */ /* for (i=statespace->zeta_modulus-1; i >= 0; i--) */ -/* { */ -/* IEz[k][i]=CreateVector(statespace->nz); */ -/* IEzz[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ -/* ISEz[k][i]=CreateVector(statespace->nz); */ -/* } */ +/* { */ +/* IEz[k][i]=CreateVector(statespace->nz); */ +/* IEzz[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ +/* ISEz[k][i]=CreateVector(statespace->nz); */ +/* } */ /* for (i=statespace->nstates-1; i >= 0; i--) */ -/* { */ -/* Ez1[k][i]=CreateVector(statespace->nz); */ -/* Ezz1[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ -/* Ez[k][i]=CreateVector(statespace->nz); */ -/* Ezz[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ -/* SEz[k][i]=CreateVector(statespace->nz); */ -/* } */ +/* { */ +/* Ez1[k][i]=CreateVector(statespace->nz); */ +/* Ezz1[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ +/* Ez[k][i]=CreateVector(statespace->nz); */ +/* Ezz[k][i]=CreateMatrix(statespace->nz,statespace->nz); */ +/* SEz[k][i]=CreateVector(statespace->nz); */ +/* } */ /* } */ /* ConditionalFilter(0,h-1,y,statespace->Ez[t0],statespace->Ezz[t0],model->V[t0],model->sv->Q, */ -/* Pxi,Pxi1,IEz,IEzz,Ez1,Ezz1,Ez,Ezz,statespace); */ +/* Pxi,Pxi1,IEz,IEzz,Ez1,Ezz1,Ez,Ezz,statespace); */ /* SmoothProbabilities_MSStateSpace(0,h-1,SPxi,Pxi,Pxi1,model->sv->Q); */ @@ -2519,34 +2519,34 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s /* { */ /* F=dw_CreateArray_vector(h); */ /* for (k=h-1; k >= 0; k--) */ -/* F[k]=CreateVector(statespace->ny); */ +/* F[k]=CreateVector(statespace->ny); */ /* } */ /* for (k=h-1; k >= 0; k--) */ /* { */ /* InitializeVector(F[k],0.0); */ /* if (statespace->zeta_modulus > statespace->nbasestates) */ -/* { */ -/* IntegrateStatesSingle(SPzeta,SPxi[k],statespace->zeta_modulus,statespace->nbasestates,2); */ -/* IntegrateStatesSingleV(z,SPzeta,ISEz[k],statespace->nbasestates,statespace->zeta_modulus/statespace->nbasestates,2); */ -/* IntegrateStatesSingle(SPs,SPzeta,statespace->nbasestates,statespace->zeta_modulus/statespace->nbasestates,2); */ -/* for (s=statespace->nbasestates-1; s >= 0; s--) */ -/* { */ -/* ProductMV(u,statespace->H[s],z[s]); */ -/* AddVV(u,statespace->a[s],u); */ -/* LinearCombinationV(F[k],1.0,F[k],ElementV(SPs,s),u); */ -/* } */ -/* } */ +/* { */ +/* IntegrateStatesSingle(SPzeta,SPxi[k],statespace->zeta_modulus,statespace->nbasestates,2); */ +/* IntegrateStatesSingleV(z,SPzeta,ISEz[k],statespace->nbasestates,statespace->zeta_modulus/statespace->nbasestates,2); */ +/* IntegrateStatesSingle(SPs,SPzeta,statespace->nbasestates,statespace->zeta_modulus/statespace->nbasestates,2); */ +/* for (s=statespace->nbasestates-1; s >= 0; s--) */ +/* { */ +/* ProductMV(u,statespace->H[s],z[s]); */ +/* AddVV(u,statespace->a[s],u); */ +/* LinearCombinationV(F[k],1.0,F[k],ElementV(SPs,s),u); */ +/* } */ +/* } */ /* else */ -/* { */ -/* IntegrateStatesSingle(SPs,SPxi[k],statespace->nbasestates,statespace->zeta_modulus,2); */ -/* for (s=statespace->nbasestates-1; s >= 0; s--) */ -/* { */ -/* ProductMV(u,statespace->H[s],ISEz[k][s]); */ -/* AddVV(u,statespace->a[s],u); */ -/* LinearCombinationV(F[k],1.0,F[k],ElementV(SPs,s),u); */ -/* } */ -/* } */ +/* { */ +/* IntegrateStatesSingle(SPs,SPxi[k],statespace->nbasestates,statespace->zeta_modulus,2); */ +/* for (s=statespace->nbasestates-1; s >= 0; s--) */ +/* { */ +/* ProductMV(u,statespace->H[s],ISEz[k][s]); */ +/* AddVV(u,statespace->a[s],u); */ +/* LinearCombinationV(F[k],1.0,F[k],ElementV(SPs,s),u); */ +/* } */ +/* } */ /* } */ /* // Clean up */ @@ -3590,19 +3590,19 @@ TMatrix* MakeA0_All(TMatrix *A0, T_VAR_Parameters *p) if (!A0) { if (!(A0=dw_CreateArray_matrix(p->nstates))) - return (TMatrix*)NULL; + return (TMatrix*)NULL; } else if ((dw_DimA(A0) != p->nstates)) { - dw_Error(SIZE_ERR); - return (TMatrix*)NULL; + dw_Error(SIZE_ERR); + return (TMatrix*)NULL; } for (s=p->nstates-1; s >= 0; s--) if (!(A0[s]=MakeA0(A0[s],s,p))) { - if (A0_in != A0) dw_FreeArray(A0); - return (TMatrix*)NULL; + if (A0_in != A0) dw_FreeArray(A0); + return (TMatrix*)NULL; } return A0; } @@ -3644,19 +3644,19 @@ TMatrix* MakeAplus_All(TMatrix *Aplus, T_VAR_Parameters *p) if (!Aplus) { if (!(Aplus=dw_CreateArray_matrix(p->nstates))) - return (TMatrix*)NULL; + return (TMatrix*)NULL; } else if ((dw_DimA(Aplus) != p->nstates)) { - dw_Error(SIZE_ERR); - return (TMatrix*)NULL; + dw_Error(SIZE_ERR); + return (TMatrix*)NULL; } for (s=p->nstates-1; s >= 0; s--) if (!(Aplus[s]=MakeAplus(Aplus[s],s,p))) { - if (Aplus_in != Aplus) dw_FreeArray(Aplus); - return (TMatrix*)NULL; + if (Aplus_in != Aplus) dw_FreeArray(Aplus); + return (TMatrix*)NULL; } return Aplus; } @@ -3694,19 +3694,19 @@ TMatrix* MakeZeta_All(TMatrix *Zeta, T_VAR_Parameters *p) if (!Zeta) { if (!(Zeta=dw_CreateArray_matrix(p->nstates))) - return (TMatrix*)NULL; + return (TMatrix*)NULL; } else if ((dw_DimA(Zeta) != p->nstates)) { - dw_Error(SIZE_ERR); - return (TMatrix*)NULL; + dw_Error(SIZE_ERR); + return (TMatrix*)NULL; } for (s=p->nstates-1; s >= 0; s--) if (!(Zeta[s]=MakeZeta(Zeta[s],s,p))) { - if (Zeta_in != Zeta) dw_FreeArray(Zeta); - return (TMatrix*)NULL; + if (Zeta_in != Zeta) dw_FreeArray(Zeta); + return (TMatrix*)NULL; } return Zeta; } diff --git a/matlab/swz/c-code/sbvar/var/VARio.c b/matlab/swz/c-code/sbvar/var/VARio.c index eb9041302..f60c370c3 100644 --- a/matlab/swz/c-code/sbvar/var/VARio.c +++ b/matlab/swz/c-code/sbvar/var/VARio.c @@ -580,36 +580,36 @@ int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model) { A0[s]=CreateMatrix(p->nvars,p->nvars); for (j=0; j < p->nvars; j++) - for (i=0; i < p->nvars; i++) - if (fscanf(f_in," %lf ",&ElementM(A0[s],i,j)) != 1) - goto ERROR; + for (i=0; i < p->nvars; i++) + if (fscanf(f_in," %lf ",&ElementM(A0[s],i,j)) != 1) + goto ERROR; Aplus[s]=CreateMatrix(p->npre,p->nvars); for (j=0; j < p->nvars; j++) - for (i=0; i < p->npre; i++) - if (fscanf(f_in," %lf ",&ElementM(Aplus[s],i,j)) != 1) - goto ERROR; + for (i=0; i < p->npre; i++) + if (fscanf(f_in," %lf ",&ElementM(Aplus[s],i,j)) != 1) + goto ERROR; Zeta[s]=CreateVector(p->nvars); for (j=0; j < p->nvars; j++) - if (fscanf(f_in," %lf ",&ElementV(Zeta[s],j)) != 1) - goto ERROR; - else - if (ElementV(Zeta[s],j) < 0.0) - goto ERROR; + if (fscanf(f_in," %lf ",&ElementV(Zeta[s],j)) != 1) + goto ERROR; + else + if (ElementV(Zeta[s],j) < 0.0) + goto ERROR; } // Set A0, Aplus, and Zeta for (j=0; j < p->nvars; j++) for (s=0; s < p->nstates; s++) { - for (i=0; i < p->nvars; i++) - ElementV(p->A0[j][p->coef_states[j][s]],i)=ElementM(A0[s],i,j); + for (i=0; i < p->nvars; i++) + ElementV(p->A0[j][p->coef_states[j][s]],i)=ElementM(A0[s],i,j); - for (i=0; i < p->npre; i++) - ElementV(p->Aplus[j][p->coef_states[j][s]],i)=ElementM(Aplus[s],i,j); + for (i=0; i < p->npre; i++) + ElementV(p->Aplus[j][p->coef_states[j][s]],i)=ElementM(Aplus[s],i,j); - p->Zeta[j][p->var_states[j][s]]=ElementV(Zeta[s],j); + p->Zeta[j][p->var_states[j][s]]=ElementV(Zeta[s],j); } // Update b0, bplus, lambda, psi diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index 0d7d90e3d..5f519e75b 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -62,49 +62,49 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste { // Read parameters and push them into model if (!posterior_file) - done=1; + done=1; else - if (!ReadBaseTransitionMatricesFlat(posterior_file,model) || !Read_VAR_ParametersFlat(posterior_file,model)) - { - done=2; - printf("total posterior draws processed - %d\n",i); - } - else - if (i++ == n) - { - printf("%d posterior draws processed\n",i); - n+=1000; - } + if (!ReadBaseTransitionMatricesFlat(posterior_file,model) || !Read_VAR_ParametersFlat(posterior_file,model)) + { + done=2; + printf("total posterior draws processed - %d\n",i); + } + else + if (i++ == n) + { + printf("%d posterior draws processed\n",i); + n+=1000; + } - if (done != 2) - { - // Get filtered probability at time T - for (j=p->nstates-1; j >= 0; j--) - ElementV(init_prob,j)=ProbabilityStateConditionalCurrent(j,T,model); + if (done != 2) + { + // Get filtered probability at time T + for (j=p->nstates-1; j >= 0; j--) + ElementV(init_prob,j)=ProbabilityStateConditionalCurrent(j,T,model); - for (k=draws; k > 0; k--) - { - // Draw time T regime - m=DrawDiscrete(init_prob); + for (k=draws; k > 0; k--) + { + // Draw time T regime + m=DrawDiscrete(init_prob); - // Draw regimes from time T+1 through T+h inclusive - for (j=0; j < h; j++) - { - ColumnVector(prob,model->sv->Q,m); - S[j]=m=DrawDiscrete(prob); - } + // Draw regimes from time T+1 through T+h inclusive + for (j=0; j < h; j++) + { + ColumnVector(prob,model->sv->Q,m); + S[j]=m=DrawDiscrete(prob); + } - // Draw shocks - for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); // InitializeVector(shocks[i],0.0); + // Draw shocks + for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); // InitializeVector(shocks[i],0.0); - // Compute forecast - if (!forecast_base(forecast,h,initial,shocks,S,model)) - goto ERROR_EXIT; + // Compute forecast + if (!forecast_base(forecast,h,initial,shocks,S,model)) + goto ERROR_EXIT; - // Accumulate impulse response - AddMatrixObservation(forecast,histogram); - } - } + // Accumulate impulse response + AddMatrixObservation(forecast,histogram); + } + } } for (i=0; i < DimV(percentiles); i++) @@ -152,7 +152,7 @@ ERROR_EXIT: */ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, - FILE *posterior_file, int s, int T, int h, TStateModel *model) + FILE *posterior_file, int s, int T, int h, TStateModel *model) { T_VAR_Parameters *p; int done=0, rtrn=0, *S, i=0, j, k, m, n=1000; @@ -186,35 +186,35 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, { // Read parameters and push them into model if (!posterior_file) - done=1; + done=1; else - if (!ReadBaseTransitionMatricesFlat(posterior_file,model) || !Read_VAR_ParametersFlat(posterior_file,model)) - { - done=2; - printf("total posterior draws processed - %d\n",i); - } - else - if (i++ == n) - { - printf("%d posterior draws processed\n",i); - n+=1000; - } + if (!ReadBaseTransitionMatricesFlat(posterior_file,model) || !Read_VAR_ParametersFlat(posterior_file,model)) + { + done=2; + printf("total posterior draws processed - %d\n",i); + } + else + if (i++ == n) + { + printf("%d posterior draws processed\n",i); + n+=1000; + } if (done != 2) - { - for (k=draws; k > 0; k--) - { - // Draw shocks - for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); // InitializeVector(shocks[i],0.0); + { + for (k=draws; k > 0; k--) + { + // Draw shocks + for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); // InitializeVector(shocks[i],0.0); - // Compute forecast - if (!forecast_base(forecast,h,initial,shocks,S,model)) - goto ERROR_EXIT; + // Compute forecast + if (!forecast_base(forecast,h,initial,shocks,S,model)) + goto ERROR_EXIT; - // Accumulate impulse response - AddMatrixObservation(forecast,histogram); - } - } + // Accumulate impulse response + AddMatrixObservation(forecast,histogram); + } + } } for (i=0; i < DimV(percentiles); i++) @@ -325,11 +325,11 @@ int main(int nargs, char **args) // specification filename if (!spec) - sprintf(spec=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + sprintf(spec=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); // parameter filename if (!parm) - sprintf(parm=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + sprintf(parm=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); } // horizon @@ -340,11 +340,11 @@ int main(int nargs, char **args) swz_fprintf_err("No specification filename given\n"); swz_fprintf_err("Command line syntax:\n" " -ft : file tag\n" - " -fs : specification filename (est_final_.dat)\n" - " -fp : parameters filename (specification filename)\n" + " -fs : specification filename (est_final_.dat)\n" + " -fp : parameters filename (specification filename)\n" " -fh : parameter header (Posterior mode: )\n" " -horizon : horizon for the forecast (12)\n" - ); + ); swzExit(1); } @@ -377,8 +377,8 @@ int main(int nargs, char **args) /* swzFree(out_filename); */ /* printf("Constructing mean forecast\n"); */ /* if (F=dw_state_space_mean_unconditional_forecast((TVector*)NULL,h,statespace->nobs,model)) */ - /* for (i=0; i < h; i++) */ - /* dw_PrintVector(f_out,F[i],"%le "); */ + /* for (i=0; i < h; i++) */ + /* dw_PrintVector(f_out,F[i],"%le "); */ /* fclose(f_out); */ /* return; */ /* } */ @@ -390,10 +390,10 @@ int main(int nargs, char **args) fmt="draws_%s.dat"; sprintf(post=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); if (!(posterior_file=fopen(post,"rt"))) - { - printf("Unable to open draws file: %s\n",post); - swzExit(0); - } + { + printf("Unable to open draws file: %s\n",post); + swzExit(0); + } // Get thinning factor from command line thin=dw_ParseInteger_String(nargs,args,"thin",1); @@ -417,66 +417,66 @@ int main(int nargs, char **args) if ((i=dw_FindArgument_String(nargs,args,"percentiles")) == -1) if (dw_FindArgument_String(nargs,args,"error_bands") == -1) { - percentiles=CreateVector(1); - ElementV(percentiles,0)=0.5; + percentiles=CreateVector(1); + ElementV(percentiles,0)=0.5; } else { - percentiles=CreateVector(3); - ElementV(percentiles,0)=0.16; ElementV(percentiles,1)=0.5; ElementV(percentiles,2)=0.84; + percentiles=CreateVector(3); + ElementV(percentiles,0)=0.16; ElementV(percentiles,1)=0.5; ElementV(percentiles,2)=0.84; } else if ((i+1 < nargs) && dw_IsInteger(args[i+1]) && ((n=atoi(args[i+1])) > 0) && (i+1+n < nargs)) { - percentiles=CreateVector(n); - for (j=0; j < n; j++) - if (!dw_IsFloat(args[i+2+j])|| ((ElementV(percentiles,j)=atof(args[i+2+j])) <= 0.0) - || (ElementV(percentiles,j) >= 1.0)) break; - if (j < n) - { - FreeVector(percentiles); - printf("forecasting command line: Error parsing percentiles\n"); - swzExit(0); - } + percentiles=CreateVector(n); + for (j=0; j < n; j++) + if (!dw_IsFloat(args[i+2+j])|| ((ElementV(percentiles,j)=atof(args[i+2+j])) <= 0.0) + || (ElementV(percentiles,j) >= 1.0)) break; + if (j < n) + { + FreeVector(percentiles); + printf("forecasting command line: Error parsing percentiles\n"); + swzExit(0); + } } else { - printf("forecasting command line(): Error parsing percentiles\n"); - swzExit(0); + printf("forecasting command line(): Error parsing percentiles\n"); + swzExit(0); } if (dw_FindArgument_String(nargs,args,"regimes") != -1) for (s=0; s < p->nstates; s++) { - rewind(posterior_file); - fmt="forecasts_percentiles_regime_%d_%s.prn"; - sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); - f_out=fopen(out_filename,"wt"); - swzFree(out_filename); - printf("Constructing percentiles for forecasts - regime %d\n",s); - forecast_percentile_regime(f_out,percentiles,draws,posterior_file,s,p->nobs,horizon,model); - fclose(f_out); + rewind(posterior_file); + fmt="forecasts_percentiles_regime_%d_%s.prn"; + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + f_out=fopen(out_filename,"wt"); + swzFree(out_filename); + printf("Constructing percentiles for forecasts - regime %d\n",s); + forecast_percentile_regime(f_out,percentiles,draws,posterior_file,s,p->nobs,horizon,model); + fclose(f_out); } else if (((s=dw_ParseInteger_String(nargs,args,"regime",-1)) >= 0) && (s < p->nstates)) { - fmt="forecasts_percentiles_regime_%d_%s.prn"; - sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); - f_out=fopen(out_filename,"wt"); - swzFree(out_filename); - printf("Constructing percentiles for forecasts - regime %d\n",s); - forecast_percentile_regime(f_out,percentiles,draws,posterior_file,s,p->nobs,horizon,model); - fclose(f_out); + fmt="forecasts_percentiles_regime_%d_%s.prn"; + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + f_out=fopen(out_filename,"wt"); + swzFree(out_filename); + printf("Constructing percentiles for forecasts - regime %d\n",s); + forecast_percentile_regime(f_out,percentiles,draws,posterior_file,s,p->nobs,horizon,model); + fclose(f_out); } else { - fmt="forecasts_percentiles_%s.prn"; - sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); - f_out=fopen(out_filename,"wt"); - swzFree(out_filename); - printf("Constructing percentiles for forecasts - %d draws of shocks/regimes per posterior value\n",draws); - forecast_percentile(f_out,percentiles,draws,posterior_file,p->nobs,horizon,model); - fclose(f_out); + fmt="forecasts_percentiles_%s.prn"; + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + f_out=fopen(out_filename,"wt"); + swzFree(out_filename); + printf("Constructing percentiles for forecasts - %d draws of shocks/regimes per posterior value\n",draws); + forecast_percentile(f_out,percentiles,draws,posterior_file,p->nobs,horizon,model); + fclose(f_out); } if (posterior_file) fclose(posterior_file); diff --git a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c index 627a18b21..7449cfda7 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -15,7 +15,7 @@ static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECIS static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size); static PRECISION Percentile(PRECISION percentile, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size); static TMatrix MakeHistogram(int low, int *h, PRECISION min, PRECISION max,int intervals, int sample_size, - PRECISION min_out, PRECISION max_out, int bins); + PRECISION min_out, PRECISION max_out, int bins); static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins); @@ -104,14 +104,14 @@ void AddMatrixObservation(TMatrix X, TMatrixHistogram *h) { for (i=h->rows-1; i >= 0; i--) for (j=h->cols-1; j >= 0; j--) - { - h->low[i][j]=h->high[i][j]=0; - for (k=h->intervals-1; k >= 0; k--) h->freq[i][j][k]=0; - } + { + h->low[i][j]=h->high[i][j]=0; + for (k=h->intervals-1; k >= 0; k--) h->freq[i][j][k]=0; + } if (h->type == HISTOGRAM_VARIABLE) for (i=h->rows-1; i >= 0; i--) - for (j=h->cols-1; j >= 0; j--) - ElementM(h->Min,i,j)=ElementM(h->Max,i,j)=ElementM(X,i,j); + for (j=h->cols-1; j >= 0; j--) + ElementM(h->Min,i,j)=ElementM(h->Max,i,j)=ElementM(X,i,j); } if (h->type == HISTOGRAM_FIXED) @@ -245,12 +245,12 @@ void AddVectorObservation(TVector x, TVectorHistogram *h) { for (i=h->dim-1; i >= 0; i--) { - h->low[i]=h->high[i]=0; - for (k=h->intervals-1; k >= 0; k--) h->freq[i][k]=0; + h->low[i]=h->high[i]=0; + for (k=h->intervals-1; k >= 0; k--) h->freq[i][k]=0; } if (h->type == HISTOGRAM_VARIABLE) for (i=h->dim-1; i >= 0; i--) - ElementV(h->Min,i)=ElementV(h->Max,i)=ElementV(x,i); + ElementV(h->Min,i)=ElementV(h->Max,i)=ElementV(x,i); } if (h->type == HISTOGRAM_FIXED) @@ -691,57 +691,57 @@ static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECI if (lo == intervals) { - min_out=max-1.0; - max_out=max+1.0; + min_out=max-1.0; + max_out=max+1.0; } else { - if (high > 0) - hi=intervals; - else - for (hi=intervals-1; !h[hi]; hi--); + if (high > 0) + hi=intervals; + else + for (hi=intervals-1; !h[hi]; hi--); - if (lo >= 0) - if (hi < intervals) - { - min_out=min+lo*inc; - max_out=min+(hi+1)*inc; - } - else - { - min_out=min+lo*inc; - if (bins == 1) - max_out=(1+SQRT_MACHINE_EPSILON)*max; - else - { - inc=(1-SQRT_MACHINE_EPSILON)*(max - min_out)/(PRECISION)(bins-1); - max_out=max + inc; - } - } - else - if (hi < intervals) - { - max_out=min+(hi+1)*inc; - if (bins == 1) - min_out=(1-SQRT_MACHINE_EPSILON)*min; - else - { - inc=(1-SQRT_MACHINE_EPSILON)*(max_out - min)/(PRECISION)(bins-1); - min_out=min - inc; - } - } - else - if (bins <= 2) - { - min_out=(1-SQRT_MACHINE_EPSILON)*min; - max_out=(1+SQRT_MACHINE_EPSILON)*max; - } - else - { - inc=(1-SQRT_MACHINE_EPSILON)*(max_out - min)/(PRECISION)(bins-2); - min_out=min - inc; - max_out=max +inc; - } + if (lo >= 0) + if (hi < intervals) + { + min_out=min+lo*inc; + max_out=min+(hi+1)*inc; + } + else + { + min_out=min+lo*inc; + if (bins == 1) + max_out=(1+SQRT_MACHINE_EPSILON)*max; + else + { + inc=(1-SQRT_MACHINE_EPSILON)*(max - min_out)/(PRECISION)(bins-1); + max_out=max + inc; + } + } + else + if (hi < intervals) + { + max_out=min+(hi+1)*inc; + if (bins == 1) + min_out=(1-SQRT_MACHINE_EPSILON)*min; + else + { + inc=(1-SQRT_MACHINE_EPSILON)*(max_out - min)/(PRECISION)(bins-1); + min_out=min - inc; + } + } + else + if (bins <= 2) + { + min_out=(1-SQRT_MACHINE_EPSILON)*min; + max_out=(1+SQRT_MACHINE_EPSILON)*max; + } + else + { + inc=(1-SQRT_MACHINE_EPSILON)*(max_out - min)/(PRECISION)(bins-2); + min_out=min - inc; + max_out=max +inc; + } } }