From 7808df09354085f7057f9591af29b16967e984b2 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 17:47:36 +0200 Subject: [PATCH 01/12] SWZ: original forecast.c, dw_histogram.c and dw_histogram.h from Dan --- matlab/swz/c-code/sbvar/var/forecast.c | 486 ++++++++++++ .../DWCcode/histogram/dw_histogram.c | 749 ++++++++++++++++++ .../DWCcode/histogram/dw_histogram.h | 77 ++ 3 files changed, 1312 insertions(+) create mode 100644 matlab/swz/c-code/sbvar/var/forecast.c create mode 100644 matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c create mode 100644 matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c new file mode 100644 index 000000000..b590e5536 --- /dev/null +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -0,0 +1,486 @@ + +#include "switch.h" +#include "switchio.h" +#include "VARio.h" +#include "dw_parse_cmd.h" +#include "dw_ascii.h" +#include "dw_histogram.h" + +#include +#include + +/* + Assumes + f_out : valid FILE pointer + percentiles : vector of numbers between 0 and 1 inclusive + draws : number of draws of shocks and regimes to make for each posterior draw + posterior_file : FILE pointer to file containing posterior draws. If null, current parameters are used. + T : last observation to treat as data. Usually equals model->nobs. + h : non-negative integer + model : point to valid TStateModel structure + + Results: + Computes and prints to the file f_out the requested percentiles for forecasts + of the observables. + + Returns: + One upon success and zero otherwise. + +*/ +int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *posterior_file, int T, int h, TStateModel *model) +{ + T_VAR_Parameters *p; + int done=0, rtrn=0, *S, i=0, j, k, m, n=1000; + TVector init_prob, prob, *shocks, initial; + TMatrix forecast; + TMatrixHistogram *histogram; + + // quick check of passed parameters + if (!f_out || !percentiles || (draws <= 0) || (T < 0) || (h < 0) || !model) return 0; + + p=(T_VAR_Parameters*)(model->theta); + + if (T > p->nobs) return 0; + + // allocate memory + S=(int*)malloc(h*sizeof(int)); + forecast=CreateMatrix(h,p->nvars); + histogram=CreateMatrixHistogram(h,p->nvars,100,HISTOGRAM_VARIABLE); + initial=CreateVector(p->npre); + shocks=dw_CreateArray_vector(h); + for (i=h-1; i >= 0; i--) shocks[i]=CreateVector(p->nvars); + init_prob=CreateVector(p->nstates); + prob=CreateVector(p->nstates); + + // Initial value + EquateVector(initial,p->X[T]); + + i=0; + while (!done) + { + // Read parameters and push them into model + if (!posterior_file) + 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 (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); + + // 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); + + // Compute forecast + if (!forecast_base(forecast,h,initial,shocks,S,model)) + goto ERROR_EXIT; + + // Accumulate impulse response + AddMatrixObservation(forecast,histogram); + } + } + } + + for (i=0; i < DimV(percentiles); i++) + { + MatrixPercentile(forecast,ElementV(percentiles,i),histogram); + dw_PrintMatrix(f_out,forecast,"%lg "); + fprintf(f_out,"\n"); + } + + rtrn=1; + +ERROR_EXIT: + FreeMatrixHistogram(histogram); + FreeMatrix(forecast); + free(S); + FreeVector(initial); + FreeVector(prob); + FreeVector(init_prob); + dw_FreeArray(shocks); + + return rtrn; +} + +/* + Assumes + f_out : valid FILE pointer + percentiles : vector of numbers between 0 and 1 inclusive + draws : number of draws of shocks to make for each posterior draw + posterior_file : FILE pointer to file containing posterior draws. If null, current parameters are used. + s : base state + T : last observation to treat as data. Usually equals model->nobs. + h : non-negative integer + model : point to valid TStateModel/T_MSStateSpace structure + + Results: + Computes and prints to the file f_out the requested percentiles for forecasts + of the observables. + + Returns: + One upon success and zero otherwise. + + Notes: + The regime at time T is drawn from the filtered probabilities at time t, and + is set to s there after. + +*/ +int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, + 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; + TVector init_prob, prob, *shocks, initial; + TMatrix forecast; + TMatrixHistogram *histogram; + + // quick check of passed parameters + if (!f_out || !percentiles || (draws <= 0) || (T < 0) || (h < 0) || !model) return 0; + + p=(T_VAR_Parameters*)(model->theta); + + if (T > p->nobs) return 0; + + // allocate memory + S=(int*)malloc(h*sizeof(int)); + for (i=0; i < h; i++) S[i]=s; + forecast=CreateMatrix(h,p->nvars); + histogram=CreateMatrixHistogram(h,p->nvars,100,HISTOGRAM_VARIABLE); + initial=CreateVector(p->npre); + shocks=dw_CreateArray_vector(h); + for (i=h-1; i >= 0; i--) shocks[i]=CreateVector(p->nvars); + init_prob=CreateVector(p->nstates); + prob=CreateVector(p->nstates); + + // Initial value + EquateVector(initial,p->X[T]); + + i=0; + while (!done) + { + // Read parameters and push them into model + if (!posterior_file) + 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 (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); + + // Compute forecast + if (!forecast_base(forecast,h,initial,shocks,S,model)) + goto ERROR_EXIT; + + // Accumulate impulse response + AddMatrixObservation(forecast,histogram); + } + } + } + + for (i=0; i < DimV(percentiles); i++) + { + MatrixPercentile(forecast,ElementV(percentiles,i),histogram); + dw_PrintMatrix(f_out,forecast,"%lg "); + fprintf(f_out,"\n"); + } + + rtrn=1; + +ERROR_EXIT: + FreeMatrixHistogram(histogram); + FreeMatrix(forecast); + free(S); + FreeVector(initial); + FreeVector(prob); + FreeVector(init_prob); + dw_FreeArray(shocks); + + return rtrn; +} + + +/* + Attempt to set up model from command line. Command line options are the + following + + -ft + If this argument exists, then the following is attempted: + specification file name = est_final_.dat + output file name = ir__regime_.dat + parameters file name = est_final_.dat + header = "Posterior mode: " + + -fs + If this argument exists, then the specification file name is . + The argument -fs takes precedence over -ft. + + -fp + If this argument exists, then the parameters file name is . The + argument -fp takes precedence over -ft. The default value is the filename + associated with the argument -fs. + + -ph
+ If this argument exists, then the header for the parameters file is +
. The default value is "Posterior mode: ". + + -horizon + If this argument exists, then the horizon of the impulse responses is given + by the passed integer. The default value is 12. + + -error_bands + Output error bands. (default = off - only median is computed) + + -percentiles n p_1 p_2 ... p_n + Percentiles to compute. The first parameter after percentiles must be the + number of percentiles and the following values are the actual percentiles. + default = 3 0.16 0.50 0.84 if error_bands flag is set + = 1 0.50 otherwise + + -parameter_uncertainty + Apply parameter uncertainty when computing error bands. + + -shocks_per_parameter + Number of shocks and regime paths to draw for each parameter draw. The + default value is 1 if parameter_uncertainty is set and 10,000 otherwise. + + -thin + Thinning factor. Only 1/thin of the draws in posterior draws file are + used. The default value is 1. + + -regimes + Produces forecasts as if each regime were permanent. (default = off) + + -regime + Produces forecasts as if regime were permanent. + + -mean + Produces mean forecast. (default = off) + +*/ +int main(int nargs, char **args) +{ + char *spec=(char*)NULL, *parm=(char*)NULL, *head=(char*)NULL, *post=(char*)NULL, *out_filename, *tag, *buffer, *fmt; + TStateModel *model; + T_VAR_Parameters *p; + TVector percentiles=(TVector)NULL; + int s, horizon, thin, draws, i, j, n; + FILE *f_out, *posterior_file; + + // specification filename + if (buffer=dw_ParseString_String(nargs,args,"fs",(char*)NULL)) + strcpy(spec=(char*)malloc(strlen(buffer)+1),buffer); + + // parameter filename + if (buffer=dw_ParseString_String(nargs,args,"fp",(char*)NULL)) + strcpy(parm=(char*)malloc(strlen(buffer)+1),buffer); + + // header + if (buffer=dw_ParseString_String(nargs,args,"ph",(char*)NULL)) + strcpy(head=(char*)malloc(strlen(buffer)+1),buffer); + + // file tag + if (tag=dw_ParseString_String(nargs,args,"ft",(char*)NULL)) + { + fmt="est_final_%s.dat"; + + // specification filename + if (!spec) + sprintf(spec=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + + // parameter filename + if (!parm) + sprintf(parm=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + } + + // horizon + horizon=dw_ParseInteger_String(nargs,args,"horizon",12); + + if (!spec) + { + fprintf(stderr,"No specification filename given\n"); + fprintf(stderr,"Command line syntax:\n" + " -ft : file tag\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" + ); + exit(1); + } + + if (!parm) + strcpy(parm=(char*)malloc(strlen(spec)+1),spec); + + if (!head) + { + buffer="Posterior mode: "; + strcpy(head=(char*)malloc(strlen(buffer)+1),buffer); + } + + model=Read_VAR_Specification((FILE*)NULL,spec); + ReadTransitionMatrices((FILE*)NULL,parm,head,model); + Read_VAR_Parameters((FILE*)NULL,parm,head,model); + p=(T_VAR_Parameters*)(model->theta); + + free(spec); + free(head); + free(parm); + + //============================= Compute forecasts ============================= + + // Mean forecast + /* if (dw_FindArgument_String(nargs,args,"mean") != -1) */ + /* { */ + /* fmt="forecasts_mean_%s.prn"; */ + /* sprintf(out_filename=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); */ + /* f_out=fopen(out_filename,"wt"); */ + /* free(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 "); */ + /* fclose(f_out); */ + /* return; */ + /* } */ + + // Parameter uncertainty + if (dw_FindArgument_String(nargs,args,"parameter_uncertainity") != -1) + { + // Open posterior draws file + fmt="draws_%s.dat"; + sprintf(post=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + if (!(posterior_file=fopen(post,"rt"))) + { + printf("Unable to open draws file: %s\n",post); + exit(0); + } + + // Get thinning factor from command line + thin=dw_ParseInteger_String(nargs,args,"thin",1); + + // Get shocks_per_parameter from command line + draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",1); + } + else + { + // Using posterior estimate + posterior_file=(FILE*)NULL; + + // thinning factor not used + thin=1; + + // Get shocks_per_parameter from command line + draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",10000); + } + + // Setup percentiles + 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; + } + else + { + 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"); + exit(0); + } + } + else + { + printf("forecasting command line(): Error parsing percentiles\n"); + exit(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*)malloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + f_out=fopen(out_filename,"wt"); + free(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*)malloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + f_out=fopen(out_filename,"wt"); + free(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*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + f_out=fopen(out_filename,"wt"); + free(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); + FreeVector(percentiles); + //============================================================================= + + return 0; +} + diff --git a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c new file mode 100644 index 000000000..6384e5948 --- /dev/null +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -0,0 +1,749 @@ + +#include +#include +#include + +#include "dw_histogram.h" +#include "dw_error.h" + +static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); +static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); +static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals); + +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); +static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins); + + +/******************************************************************************* + The following set of routines create a matrix of histograms on the fly. +*******************************************************************************/ +/* + Assumes + rows > 0 + cols > 0 + intervals > 0 + type = HISTOGRAM_FIXED or HISTOGRAM_VARIABLE + + Results + Creates and returns a matrix histogram data structure. The size of the + matrix is m x n and the number of intervals is intrvls. + +*/ +TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type) +{ + int i, j; + TMatrixHistogram *h; + + if (!(h=(TMatrixHistogram *)malloc(sizeof(TMatrixHistogram)))) dw_Error(MEM_ERR); + + if (!(h->freq=(int***)malloc(rows*sizeof(int**)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + { + if (!(h->freq[i]=(int**)malloc(cols*sizeof(int*)))) dw_Error(MEM_ERR); + for (j=cols-1; j >= 0; j--) + if (!(h->freq[i][j]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + } + + if (!(h->low=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + if (!(h->low[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + + if (!(h->high=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + + h->Min=CreateMatrix(rows,cols); + h->Max=CreateMatrix(rows,cols); + + h->rows=rows; + h->cols=cols; + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h) +{ + EquateMatrix(h->Min,Min); + EquateMatrix(h->Max,Max); + h->sample_size=0; +} + +void FreeMatrixHistogram(TMatrixHistogram *h) +{ + int i, j; + for (i=h->rows-1; i >= 0; i--) + { + for (j=h->cols-1; j >= 0; j--) free(h->freq[i][j]); + free(h->freq[i]); + } + free(h->freq); + for (i=h->rows-1; i >= 0; i--) free(h->low[i]); + free(h->low); + for (i=h->rows-1; i >= 0; i--) free(h->high[i]); + free(h->high); + FreeMatrix(h->Min); + FreeMatrix(h->Max); + free(h); +} + +void AddMatrixObservation(TMatrix X, TMatrixHistogram *h) +{ + int i, j, k; + + if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); + + if (h->sample_size <= 0) + { + 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; + } + 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); + } + + if (h->type == HISTOGRAM_FIXED) + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + AddObservationFixed(ElementM(X,i,j),h->low[i]+j,h->freq[i][j],h->high[i]+j,ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals); + else + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + AddObservationVariable(ElementM(X,i,j),h->freq[i][j],&ElementM(h->Min,i,j),&ElementM(h->Max,i,j),h->intervals); + + h->sample_size++; +} + +void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h) +{ + int i, j; + + if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); + + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + ElementM(X,i,j)=Percentile(percentile,h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + For 0 <= i < h->rows and 0 <= j < h->cols, let + + I[i][j][k]=(h->min[i][j] + k*inc[i][j], h->min[i][j] + (k+1)*inc[i][j]), + + where inc[i][j]=(h->max[i][j] - h->min[i][j])/h->samples_size. The + distribution is uniform on I[i][k][j] and + + P(h->min[i][j] + k*inc[i][j] < x[i][j] < h->min[i][j] + (k+1)*inc[i][j]) + = h->freq[i][j][k]/h->sample_size. + + Furthermore, + + P(x[i][j] < h->min[i][j]) = 0 and P(x[i][j] > h->min[i][j]) = 0. + + In addition, if h->type == FIXED, then + + P(x[i][j] = h->min[i][j]) = h->low[i][j]/h->sample_size + + and + + P(x[i][j] = h->min[i][j]) = h->high[i][j]/h->sample_size. +*/ +void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h) +{ + int i, j; + + if ((h->rows != RowM(P)) || (h->cols != ColM(P)) || + (h->rows != RowM(Level)) || (h->cols != ColM(Level))) + dw_Error(SIZE_ERR); + + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + ElementM(P,i,j)=Cumulative(ElementM(Level,i,j),h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); +} + +TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h) +{ + return MakeHistogramAuto(h->low[i][j],h->freq[i][j],h->high[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,bins); +} + +TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h) +{ + return MakeHistogram(h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,min,max,bins); +} + +/******************************************************************************* + The following set of routines create a vector of histograms on the fly. +*******************************************************************************/ +TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type) +{ + int i; + TVectorHistogram *h; + + if (!(h=(TVectorHistogram *)malloc(sizeof(TVectorHistogram)))) + dw_Error(MEM_ERR); + + if (!(h->freq=(int**)malloc(dim*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=dim-1; i >= 0; i--) + if (!(h->freq[i]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + + if (!(h->low=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->high=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + + h->Min=CreateVector(dim); + h->Max=CreateVector(dim); + + h->dim=dim; + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h) +{ + EquateVector(h->Min,Min); + EquateVector(h->Max,Max); + h->sample_size=0; +} + +void FreeVectorHistogram(TVectorHistogram *h) +{ + int i; + for (i=h->dim-1; i >= 0; i--) free(h->freq[i]); + free(h->freq); + free(h->low); + free(h->high); + FreeVector(h->Min); + FreeVector(h->Max); + free(h); +} + +void AddVectorObservation(TVector x, TVectorHistogram *h) +{ + int i, k; + + if (h->dim != DimV(x)) dw_Error(SIZE_ERR); + + if (h->sample_size <= 0) + { + 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; + } + if (h->type == HISTOGRAM_VARIABLE) + for (i=h->dim-1; i >= 0; i--) + ElementV(h->Min,i)=ElementV(h->Max,i)=ElementV(x,i); + } + + if (h->type == HISTOGRAM_FIXED) + for (i=h->dim-1; i >= 0; i--) + AddObservationFixed(ElementV(x,i),h->low+i,h->freq[i],h->high+i,ElementV(h->Min,i),ElementV(h->Max,i),h->intervals); + else + for (i=h->dim-1; i >= 0; i--) + AddObservationVariable(ElementV(x,i),h->freq[i],&ElementV(h->Min,i),&ElementV(h->Max,i),h->intervals); + + h->sample_size++; +} + +void VectorPercentile(TVector x, PRECISION percentile, TVectorHistogram *h) +{ + int i; + + if (h->dim != DimV(x)) dw_Error(SIZE_ERR); + + for (i=h->dim-1; i >= 0; i--) + ElementV(x,i)=Percentile(percentile,h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); +} + + + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + For 0 <= i < h->dim, let + + I[i][k]=(h->min[i] + k*inc[i], h->min[i] + (k+1)*inc[i]), + + where inc[i]=(h->max[i] - h->min[i])/h->samples_size. The distribution + is uniform on I[i][k] and + + P(h->min[i] + k*inc[i] < x[i] < h->min[i] + (k+1)*inc[i]) + = h->freq[i][k]/h->sample_size. + + Furthermore, + + P(x[i] < h->min[i]) = 0 and P(x[i] > h->min[i]) = 0. + + In addition, if h->type == FIXED, then + + P(x[i] = h->min[i]) = h->low[i]/h->sample_size + + and + + P(x[i] = h->min[i]) = h->high[i]/h->sample_size. +*/ +void VectorCumulative(TVector p, TVector level, TVectorHistogram *h) +{ + int i; + + if (h->dim != DimV(p) || (h->dim != DimV(level))) + dw_Error(SIZE_ERR); + + for (i=h->dim-1; i >= 0; i--) + ElementV(p,i)=Cumulative(ElementV(level,i),h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); +} + +TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h) +{ + return MakeHistogramAuto(h->low[i],h->freq[i],h->high[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,bins); +} + +TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h) +{ + return MakeHistogram(h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,min,max,bins); +} +/******************************************************************************* + The following set of routines create a scalar histogram on the fly. +*******************************************************************************/ +/* + Assumes + + Results + Creates and returns a scalar histogram data structure. +*/ +TScalarHistogram *CreateScalarHistogram(int intervals, int type) +{ + TScalarHistogram *h; + + if (!(h=(TScalarHistogram *)malloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR); + + if (!(h->freq=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h) +{ + h->Min=Min; + h->Max=Max; + h->sample_size=0; +} + +void FreeScalarHistogram(TScalarHistogram *h) +{ + free(h->freq); + free(h); +} + +void AddScalarObservation(PRECISION x, TScalarHistogram *h) +{ + int k; + + if (h->sample_size <= 0) + { + h->low=h->high=0; + for (k=h->intervals-1; k >= 0; k--) h->freq[k]=0; + if (h->type == HISTOGRAM_VARIABLE) h->Min=h->Max=x; + } + + if (h->type == HISTOGRAM_FIXED) + AddObservationFixed(x,&(h->low),h->freq,&(h->high),h->Min,h->Max,h->intervals); + else + AddObservationVariable(x,h->freq,&(h->Min),&(h->Max),h->intervals); + + h->sample_size++; +} + +PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h) +{ + return Percentile(percentile,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + Let + + I[k]=(h->min + k*inc, h->min + (k+1)*inc), + + where inc=(h->max - h->min)/h->samples_size. The distribution + is uniform on I[k] and + + P(h->min + k*inc < x < h->min + (k+1)*inc) = h->freq[k]/h->sample_size. + + Furthermore, + + P(x < h->min) = 0 and P(x > h->min) = 0. + + In addition, if h->type == FIXED, then + + P(x = h->min) = h->low/h->sample_size + + and + + P(x = h->min) = h->high/h->sample_size. +*/ +PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h) +{ + return Cumulative(level,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); +} + +TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h) +{ + return MakeHistogramAuto(h->low,h->freq,h->Min,h->high,h->Max,h->intervals,h->sample_size,bins); +} + +TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h) +{ + return MakeHistogram(h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size,min,max,bins); +} + +/*******************************************************************************/ +/***************************** Low Level Routines ******************************/ +/*******************************************************************************/ +/* + Resizes the histogram. After resizing, it is guaranteed that *min <= x <= *max. + The type of the histogram must be HISTOGRAM_VARIABLE. +*/ +static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) +{ + int i, j, k, m; + if (x > *max) + if (x - *min >= (PRECISION)intervals*(*max - *min)) + { + for (i=1; i < intervals; i++) + { + h[0]+=h[i]; + h[i]=0; + } + *max=x; + } + else + { + m=(int)ceil((x - *min)/(*max - *min)); + for (i=j=0; i < intervals; j++) + for(h[j]=h[i++], k=1; (k < m) && (i < intervals); k++) + h[j]+=h[i++]; + for ( ; j < intervals; j++) h[j]=0; + *max=*min + m*(*max - *min); + if (x > *max) *max=x; + } + else + if (x < *min) + if (*max - x >= (PRECISION)intervals*(*max - *min)) + { + for (j=intervals-1, i=intervals-2; i >= 0; i--) + { + h[j]+=h[i]; + h[i]=0; + } + *min=x; + } + else + { + m=(int)ceil((*max - x)/(*max - *min)); + for (i=j=intervals-1; i >= 0; j--) + for(h[j]=h[i--], k=1; (k < m) && (i >= 0); k++) + h[j]+=h[i--]; + for ( ; j >= 0; j--) h[j]=0; + *min=*max - m*(*max - *min); + if (x < *min) *min=x; + } +} + +/* + Adds a observation to the histogram. The type of the histogram must + be HISTOGRAM_VARIABLE. +*/ +static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) +{ + int i; + + if ((x < *min) || (x > *max)) Resize(x,h,min,max,intervals); + + if (*max > *min) + { + i=(int)(intervals*(x - *min)/(*max - *min)); + h[(i < intervals) ? i : intervals-1]++; + } + else + h[0]++; +} + +/* + Adds a observation to the histogram. The type of the histogram must + be HISTOGRAM_FIXED. +*/ +static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals) +{ + PRECISION y=floor(intervals*(x - min)/(max - min)); + if (y < 0) + (*low)++; + else + if (y < intervals) + h[(int)y]++; + else + (*high)++; +} + +/******************************************************************************/ +/******************************************************************************/ +/******************************************************************************/ + +/* + Returns the level such that the probability of observing an observation + less than or equal to level is percentile. If there is a point mass at + x, and P(y < x) <= percentile <= P(y <= x), then x is returned. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point. +*/ +static PRECISION Percentile(PRECISION percentile, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) +{ + int i; + percentile=percentile*sample_size - low; + if (percentile <= 0) return min; + for (i=0; i < intervals; i++) + if (h[i] && (percentile-=h[i]) <= 0) + return min + ((PRECISION)(i+1) + percentile/(PRECISION)h[i])*(max - min)/(PRECISION)intervals; + return max; +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also, if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point +*/ +static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) +{ + PRECISION inc=(max-min)/(PRECISION)intervals; + int i, count; + + if (level < min) return 0.0; + if (level >= max) return 1.0; + + for (count=low, i=0; i < intervals; count+=h[i++]) + if ((min+=inc) >= level) + return ((PRECISION)count + (PRECISION)h[i]*(level - min + inc)/inc)/(PRECISION)sample_size; + return 1.0; +} + +/* + Returns a histogram over the interval I=[min_out,max_out]. The matrix returned + has bins rows and 2 columns. If inc=(max_out - min_out)/bins, then the first + element of the ith row is + + min + (i + 0.5)*inc, + + which is the mid-point of the ith interval. The second element is + + P(min + i*inc < x <= min + (i + 1)*inc)/inc, + + which is the average density over the ith interval. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point. +*/ +static TMatrix MakeHistogram(int low, int *h, PRECISION min, PRECISION max, + int intervals, int sample_size, PRECISION min_out, PRECISION max_out, int bins) +{ + int i; + PRECISION inc, x, cdf_lower, cdf_upper; + TMatrix X; + + inc=(max_out-min_out)/(PRECISION)bins; + + if (inc > 0) + { + X=CreateMatrix(bins,2); + x=min_out+inc; + + cdf_lower=Cumulative(min_out,low,h,min,max,intervals,sample_size); + + for (i=0; i < bins; i++) + { + cdf_upper=Cumulative(x,low,h,min,max,intervals,sample_size); + + ElementM(X,i,0)=x - 0.5*inc; + ElementM(X,i,1)=(cdf_upper-cdf_lower)/inc; + + cdf_lower=cdf_upper; + x+=inc; + } + } + else + return (TMatrix)NULL; + + return X; +} + +/* + Automatically chooses lenth of interval over which to produce histogram and + then calls MakeHistogram(). +*/ +static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins) +{ + PRECISION inc=(max-min)/intervals, max_out, min_out; + int lo, hi; + + if ((low == sample_size) || (inc <= 0)) + { + min_out=min-1.0; + max_out=min+1.0; + } + else + { + if (low > 0) + lo=-1; + else + for (lo=0; (lo < intervals) && !h[lo]; lo++); + + if (lo == intervals) + { + 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 (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; + } + } + } + + return MakeHistogram(low,h,min,max,intervals,sample_size,min_out,max_out,bins); +} + + diff --git a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h new file mode 100644 index 000000000..801a8859d --- /dev/null +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h @@ -0,0 +1,77 @@ +#ifndef __HISTOGRAMS__ +#define __HISTOGRAMS__ + +#include "matrix.h" + +#define HISTOGRAM_FIXED 1 +#define HISTOGRAM_VARIABLE 2 + +/* Matrix histograms */ +typedef struct +{ + TMatrix Min; + TMatrix Max; + int **low; + int **high; + int ***freq; + int rows; + int cols; + int intervals; + int sample_size; + int type; +} TMatrixHistogram; + +/* Vector histograms */ +typedef struct +{ + TVector Min; + TVector Max; + int *low; + int *high; + int **freq; + int dim; + int intervals; + int sample_size; + int type; +} TVectorHistogram; + +/* Scalar histograms */ +typedef struct +{ + PRECISION Min; + PRECISION Max; + int low; + int high; + int *freq; + int intervals; + int sample_size; + int type; +} TScalarHistogram; + +TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type); +void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h); +void FreeMatrixHistogram(TMatrixHistogram *h); +void AddMatrixObservation(TMatrix X, TMatrixHistogram *h); +void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h); +void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h); +TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h); +TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h); + +TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type); +void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h); +void FreeVectorHistogram(TVectorHistogram *h); +void AddVectorObservation(TVector X, TVectorHistogram *h); +void VectorPercentile(TVector X, PRECISION percentile, TVectorHistogram *h); +void VectorCumulative(TVector p, TVector level, TVectorHistogram *h); +TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h); +TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h); + +TScalarHistogram *CreateScalarHistogram(int intervals, int type); +void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h); +void FreeScalarHistogram(TScalarHistogram *h); +void AddScalarObservation(PRECISION x, TScalarHistogram *h); +PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h); +PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h); +TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h); +TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h); +#endif From 248cc388afb2a169c4e13b8625b5b332fb546493 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 1 Sep 2010 10:43:25 +0200 Subject: [PATCH 02/12] SWZ: diffed and patched updated files from Dan --- matlab/swz/c-code/sbvar/switching/switchio.c | 77 +++++ matlab/swz/c-code/sbvar/switching/switchio.h | 2 + matlab/swz/c-code/sbvar/var/VARbase.c | 286 +++++++++++++++++++ matlab/swz/c-code/sbvar/var/VARbase.h | 9 + matlab/swz/c-code/sbvar/var/VARio.c | 68 +++++ matlab/swz/c-code/sbvar/var/VARio.h | 1 + 6 files changed, 443 insertions(+) diff --git a/matlab/swz/c-code/sbvar/switching/switchio.c b/matlab/swz/c-code/sbvar/switching/switchio.c index 90e9975d8..179b8225e 100644 --- a/matlab/swz/c-code/sbvar/switching/switchio.c +++ b/matlab/swz/c-code/sbvar/switching/switchio.c @@ -561,6 +561,9 @@ int ReadBaseTransitionMatrices_SV(FILE *f_in, TMarkovStateVariable* sv, char *he /* // Convert base transition matrix to full transition matrix. ansi-c*/ ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded); + /* Free Q */ + FreeMatrix(Q); + /* // Update ansi-c*/ if (!Update_B_from_Q_SV(sv)) { @@ -570,6 +573,7 @@ int ReadBaseTransitionMatrices_SV(FILE *f_in, TMarkovStateVariable* sv, char *he return 1; } + FreeMatrix(Q); return 0; } } @@ -712,6 +716,74 @@ void WriteBaseTransitionMatricesFlat_Headers_SV(FILE *f_out, TMarkovStateVariabl } } +/* + Returns 1 upon success and 0 upon failure. +*/ +int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) +{ + int i, j; + TMatrix Q; + PRECISION sum; + + 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; + } + 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; + } + + // 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; + } + + // Convert base transition matrix to full transition matrix. + ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded); + + // Free Q + FreeMatrix(Q); + + // Update + if (!Update_B_from_Q_SV(sv)) + { + dw_UserError("Transition matrices do not satisfy restrictions"); + return 0; + } + + return 1; + } + + return 1; +} + /* Returns 1 upon success and 0 upon failure. @@ -894,6 +966,11 @@ int WriteBaseTransitionMatrices(FILE *f, char *filename, char *header, TStateMod return rtrn; } +int ReadBaseTransitionMatricesFlat(FILE *f, TStateModel *model) +{ + return f ? ReadBaseTransitionMatricesFlat_SV(f,model->sv) : 0; +} + int WriteBaseTransitionMatricesFlat(FILE *f, TStateModel *model, char *fmt) { return f ? WriteBaseTransitionMatricesFlat_SV(f,model->sv,fmt) : 0; diff --git a/matlab/swz/c-code/sbvar/switching/switchio.h b/matlab/swz/c-code/sbvar/switching/switchio.h index c5bdbd4d8..08f7d3bcc 100644 --- a/matlab/swz/c-code/sbvar/switching/switchio.h +++ b/matlab/swz/c-code/sbvar/switching/switchio.h @@ -12,6 +12,7 @@ int WriteTransitionMatrices_SV(FILE *f_out, TMarkovStateVariable* sv, char *head int ReadBaseTransitionMatrices_SV(FILE *f_out, TMarkovStateVariable *sv, char *header, char *idstring); int WriteBaseTransitionMatrices_SV(FILE *f_out, TMarkovStateVariable *sv, char *header, char *idstring); +int ReadBaseTransitionMatricesFlat_SV(FILE *f_out, TMarkovStateVariable *sv); int WriteBaseTransitionMatricesFlat_SV(FILE *f_out, TMarkovStateVariable *sv, char *fmt); void WriteBaseTransitionMatricesFlat_Headers_SV(FILE *f_out, TMarkovStateVariable* sv, char *idstring); @@ -28,6 +29,7 @@ int WriteStates(FILE *f, char *filename, char *header, TStateModel *model); int ReadBaseTransitionMatrices(FILE *f, char *filename, char *header, TStateModel *model); int WriteBaseTransitionMatrices(FILE *f, char *filename, char *header, TStateModel *model); +int ReadBaseTransitionMatricesFlat(FILE *f, TStateModel *model); int WriteBaseTransitionMatricesFlat(FILE *f, TStateModel *model, char *fmt); /* diff --git a/matlab/swz/c-code/sbvar/var/VARbase.c b/matlab/swz/c-code/sbvar/var/VARbase.c index 21d1f7961..ad776c50a 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.c +++ b/matlab/swz/c-code/sbvar/var/VARbase.c @@ -2389,6 +2389,203 @@ void Update_lambda_psi_from_bplus(T_VAR_Parameters *p) /******************************************************************************/ /******************************************************************************/ + +/******************************************************************************/ +/********************************* Forecasts **********************************/ +/******************************************************************************/ +/* + Assumes: + forecast : horizon x nvars matrix or null pointer + horizon : positive integer - forecast horizon + initial : initial value of predetermined variables. + shocks : array of length horizon of shocks or null pointer. If null + pointer, then the shocks are all zero. Each vector is of length + nvar. + S : array of length horizon. S[t] is the state at time T+1+t. + model : pointer to valid TStateModel structure. + + Results: + Computes forecast + + Returns: + The matrix forecast upon success and null upon failure. If forecast is + null, then it created. +*/ +TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *shocks, int *S, TStateModel *model) +{ + T_VAR_Parameters *p=(T_VAR_Parameters*)(model->theta); + TMatrix *A0, *Aplus; + TVector x, y; + int i, t; + + // allocate forecast if necessary + if (!forecast && !(forecast=CreateMatrix(horizon,p->nvars))) + return (TMatrix)NULL; + + // allocate memory + y=CreateVector(p->nvars); + x=EquateVector((TVector)NULL,initial); + A0=MakeA0_All((TMatrix*)NULL,p); + Aplus=MakeAplus_All((TMatrix*)NULL,p); + + // forecast + for (t=0; t < horizon; t++) + { + 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]]]); + } + ProductInverseVM(y,y,A0[S[t]]); + for (i=p->nvars-1; i >= 0; 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)); + } + + // free memory + dw_FreeArray(Aplus); + dw_FreeArray(A0); + FreeVector(x); + FreeVector(y); + + return forecast; +} + +/* + For 1 <= k < h, y[k][i] is null if the ith coordinate of y(t0+1+k) is + unrestricted and is its value otherwise. In general, t0 is the last index for + which we have full information. It must be the case that t0 <= nobs. +*/ +/* TVector* dw_state_space_mean_conditional_forecast(TVector *F, PRECISION ***y, int h, int t0, TStateModel *model) */ +/* { */ +/* T_MSStateSpace *statespace=(T_MSStateSpace*)(model->theta); */ +/* TVector Pxi_i, *Pxi, *Pxi1, *SPxi, *Ez_i, **IEz, **Ez1, **Ez, **SEz, **ISEz, SPzeta, SPs, u, *z; */ +/* TMatrix Q, *Ezz_i, **IEzz, **Ezz1, **Ezz; */ +/* int i, k, s; */ + +/* if ((t0 > statespace->t0) && !Filter(t0,model)) return (TVector*)NULL; */ +/* if ((t0 > model->t0) && !ForwardRecursion(t0,model)) return (TVector*)NULL; */ +/* if (((t0 < model->sv->t0) || (model->sv->t1 < t0)) && !sv_ComputeTransitionMatrix(t0,model->sv,model)) return (TVector*)NULL; */ + +/* Pxi=dw_CreateArray_vector(h); */ +/* Pxi1=dw_CreateArray_vector(h); */ +/* SPxi=dw_CreateArray_vector(h); */ +/* IEz=dw_CreateRectangularArray_vector(h,statespace->zeta_modulus); */ +/* IEzz=dw_CreateRectangularArray_matrix(h,statespace->zeta_modulus); */ +/* Ez1=dw_CreateRectangularArray_vector(h,statespace->nstates); */ +/* Ezz1=dw_CreateRectangularArray_matrix(h,statespace->nstates); */ +/* Ez=dw_CreateRectangularArray_vector(h,statespace->nstates); */ +/* Ezz=dw_CreateRectangularArray_matrix(h,statespace->nstates); */ +/* SEz=dw_CreateRectangularArray_vector(h,statespace->nstates); */ +/* ISEz=dw_CreateRectangularArray_vector(h,statespace->zeta_modulus); */ +/* for (k=h-1; k >= 0; k--) */ +/* { */ +/* Pxi[k]=CreateVector(statespace->nstates); */ +/* 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); */ +/* } */ +/* 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); */ +/* } */ +/* } */ + +/* 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); */ + +/* SmoothProbabilities_MSStateSpace(0,h-1,SPxi,Pxi,Pxi1,model->sv->Q); */ + +/* SmoothMean_MSStateSpace(0,h-1,SEz,ISEz,Ez1,Ezz1,IEz,IEzz,SPxi,statespace); */ + +/* SPzeta=CreateVector(statespace->zeta_modulus); */ +/* SPs=CreateVector(statespace->nbasestates); */ +/* u=CreateVector(statespace->ny); */ +/* z=dw_CreateArray_vector(statespace->nbasestates); */ +/* for (s=statespace->nbasestates-1; s >= 0; s--) z[s]=CreateVector(statespace->nz); */ + +/* if (!F) */ +/* { */ +/* F=dw_CreateArray_vector(h); */ +/* for (k=h-1; k >= 0; k--) */ +/* 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); */ +/* } */ +/* } */ +/* 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); */ +/* } */ +/* } */ +/* } */ + +/* // Clean up */ +/* dw_FreeArray(z); */ +/* FreeVector(u); */ +/* FreeVector(SPs); */ +/* FreeVector(SPzeta); */ +/* dw_FreeArray(ISEz); */ +/* dw_FreeArray(SEz); */ +/* dw_FreeArray(Ezz); */ +/* dw_FreeArray(Ez); */ +/* dw_FreeArray(Ezz1); */ +/* dw_FreeArray(Ez1); */ +/* dw_FreeArray(IEzz); */ +/* dw_FreeArray(IEz); */ +/* dw_FreeArray(SPxi); */ +/* dw_FreeArray(Pxi1); */ +/* dw_FreeArray(Pxi); */ + +/* return F; */ +/* } */ + +/* + +*/ +/* TVector* dw_state_space_mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model) */ +/* { */ +/* T_MSStateSpace *statespace=(T_MSStateSpace*)(model->theta); */ +/* PRECISION ***y=(PRECISION***)dw_CreateMultidimensionalArrayList_scalar(3,h,statespace->ny,1); */ +/* int i, j; */ +/* for (i=h-1; i >= 0; i--) */ +/* for (j=statespace->ny-1; j >= 0; j--) */ +/* { dw_FreeArray(y[i][j]); y[i][j]=(PRECISION*)NULL; } */ +/* F=dw_state_space_mean_conditional_forecast(F,y,h,t0,model); */ +/* dw_FreeArray(y); */ +/* return F; */ +/* } */ + + /******************************************************************************/ /** Impulse Response Routines **/ /******************************************************************************/ @@ -3381,6 +3578,36 @@ TMatrix MakeA0(TMatrix A0, int s, T_VAR_Parameters *p) return A0; } +/* + Assumes: + A0 : Matrix array of length n_states or null pointer. A0[s] is either + p->nvars x p->nvars matrix or null pointer +*/ +TMatrix* MakeA0_All(TMatrix *A0, T_VAR_Parameters *p) +{ + int s; + TMatrix *A0_in=A0; + if (!A0) + { + if (!(A0=dw_CreateArray_matrix(p->nstates))) + return (TMatrix*)NULL; + } + else + if ((dw_DimA(A0) != p->nstates)) + { + 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; + } + return A0; +} + + /* Assumes: Aplus : p->npre x p->nvars matrix or null pointer @@ -3405,6 +3632,36 @@ TMatrix MakeAplus(TMatrix Aplus, int k, T_VAR_Parameters *p) return Aplus; } +/* + Assumes: + Aplus : Matrix array of length n_states or null pointer. Aplus[s] is either + p->npre x p->nvars matrix or null pointer +*/ +TMatrix* MakeAplus_All(TMatrix *Aplus, T_VAR_Parameters *p) +{ + int s; + TMatrix *Aplus_in=Aplus; + if (!Aplus) + { + if (!(Aplus=dw_CreateArray_matrix(p->nstates))) + return (TMatrix*)NULL; + } + else + if ((dw_DimA(Aplus) != p->nstates)) + { + 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; + } + return Aplus; +} + + TMatrix MakeZeta(TMatrix Zeta, int k, T_VAR_Parameters *p) { int j; @@ -3425,6 +3682,35 @@ TMatrix MakeZeta(TMatrix Zeta, int k, T_VAR_Parameters *p) return Zeta; } +/* + Assumes: + Zeta : Matrix array of length n_states or null pointer. Zeta[s] is either + p->vars x p->nvars matrix or null pointer +*/ +TMatrix* MakeZeta_All(TMatrix *Zeta, T_VAR_Parameters *p) +{ + int s; + TMatrix *Zeta_in=Zeta; + if (!Zeta) + { + if (!(Zeta=dw_CreateArray_matrix(p->nstates))) + return (TMatrix*)NULL; + } + else + if ((dw_DimA(Zeta) != p->nstates)) + { + 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; + } + return Zeta; +} + /* Assumes X: n x n matrix or null pointer in column major format diff --git a/matlab/swz/c-code/sbvar/var/VARbase.h b/matlab/swz/c-code/sbvar/var/VARbase.h index 6e8e815a2..e221daa1c 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.h +++ b/matlab/swz/c-code/sbvar/var/VARbase.h @@ -198,14 +198,23 @@ void DrawAplus(TStateModel *model); void Draw_psi(TStateModel *model); void Draw_lambda(TStateModel *model); +/* Forecasts */ +TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *shocks, int *S, TStateModel *model); + +//TVector* mean_conditional_forecast(TVector *F, PRECISION ***y, int h, int t0, TStateModel *model); +//TVector* mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model); + /* Utilities */ void ComputeDotProducts_All(T_VAR_Parameters *p); void ComputeLogAbsDetA0_All(T_VAR_Parameters *p); void ComputeLogAbsDetA0(int j, int k, T_VAR_Parameters *p); TMatrix MakeA0(TMatrix A0, int k, T_VAR_Parameters *p); +TMatrix* MakeA0_All(TMatrix *A0, T_VAR_Parameters *p); TMatrix MakeAplus(TMatrix Aplus, int k, T_VAR_Parameters *p); +TMatrix* MakeAplus_All(TMatrix *Aplus, T_VAR_Parameters *p); TMatrix MakeZeta(TMatrix Zeta, int k, T_VAR_Parameters *p); +TMatrix* MakeZeta_All(TMatrix *Zeta, T_VAR_Parameters *p); TMatrix ConstructMatrixFromColumns(TMatrix X, TVector **, int k); void UpdateStateDependentFields(T_VAR_Parameters *p, int *S); diff --git a/matlab/swz/c-code/sbvar/var/VARio.c b/matlab/swz/c-code/sbvar/var/VARio.c index 3e000128b..eb9041302 100644 --- a/matlab/swz/c-code/sbvar/var/VARio.c +++ b/matlab/swz/c-code/sbvar/var/VARio.c @@ -563,6 +563,74 @@ int Write_VAR_ParametersFlat_Headers(FILE *f_out, TStateModel *model) return 1; } +int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model) +{ + TMatrix *A0, *Aplus; + TVector *Zeta; + int i, j, s, rtrn=0; + T_VAR_Parameters *p=(T_VAR_Parameters*)(model->theta); + + // Allocate memory + A0=dw_CreateArray_matrix(p->nstates); + Aplus=dw_CreateArray_matrix(p->nstates); + Zeta=dw_CreateArray_vector(p->nstates); + + // Read File + for (s=0; s < p->nstates; s++) + { + 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; + + 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; + + 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; + } + + // 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->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); + } + + // Update b0, bplus, lambda, psi + Update_b0_bplus_from_A0_Aplus(p); + if ((p->Specification & SPEC_SIMS_ZHA) == SPEC_SIMS_ZHA) Update_lambda_psi_from_bplus(p); + + // Flags and notification that the VAR parameters have changed + p->valid_parameters=1; + ThetaChanged(model); + rtrn=1; + + ERROR: + + // Free memory + dw_FreeArray(A0); + dw_FreeArray(Aplus); + dw_FreeArray(Zeta); + + return rtrn; +} + /* For each state the VAR parameters are printed as follows A0 (by columns) diff --git a/matlab/swz/c-code/sbvar/var/VARio.h b/matlab/swz/c-code/sbvar/var/VARio.h index 942fa053c..2805beecc 100644 --- a/matlab/swz/c-code/sbvar/var/VARio.h +++ b/matlab/swz/c-code/sbvar/var/VARio.h @@ -10,6 +10,7 @@ TStateModel* Read_VAR_Specification(FILE *f, char *filename); int Write_VAR_Parameters(FILE *f, char *filename, char *id, TStateModel *model); int Read_VAR_Parameters(FILE *f, char *filename, char *id, TStateModel *model); +int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model); int Write_VAR_ParametersFlat(FILE *f, TStateModel *model, char *fmt); int Write_VAR_ParametersFlat_Headers(FILE *f_out, TStateModel *model); int Write_VAR_ParametersFlat_A0_Diagonal_One(FILE *f, TStateModel *model, char *fmt); From a5ed818758ce33b6ec5ad69ba71d0f4f7e3b35a6 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 17:55:53 +0200 Subject: [PATCH 03/12] SWZ: remove extra whitespace --- matlab/swz/c-code/sbvar/switching/switchio.c | 6 +-- matlab/swz/c-code/sbvar/var/VARbase.c | 16 +++--- matlab/swz/c-code/sbvar/var/forecast.c | 50 +++++++++---------- .../DWCcode/histogram/dw_histogram.c | 38 +++++++------- 4 files changed, 55 insertions(+), 55 deletions(-) diff --git a/matlab/swz/c-code/sbvar/switching/switchio.c b/matlab/swz/c-code/sbvar/switching/switchio.c index 179b8225e..3f58083eb 100644 --- a/matlab/swz/c-code/sbvar/switching/switchio.c +++ b/matlab/swz/c-code/sbvar/switching/switchio.c @@ -742,11 +742,11 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) 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--) + for (sum=0.0, i=sv->nbasestates-1; i >= 0; i--) if (ElementM(Q,i,j) < 0.0) { FreeMatrix(Q); @@ -761,7 +761,7 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) dw_UserError("Transition matrix columns must sum to one."); return 0; } - for (sum=1.0/sum, i=sv->nbasestates-1; i >= 0; i--) + for (sum=1.0/sum, i=sv->nbasestates-1; i >= 0; i--) ElementM(Q,i,j)*=sum; } diff --git a/matlab/swz/c-code/sbvar/var/VARbase.c b/matlab/swz/c-code/sbvar/var/VARbase.c index ad776c50a..bbed86e84 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.c +++ b/matlab/swz/c-code/sbvar/var/VARbase.c @@ -2398,8 +2398,8 @@ void Update_lambda_psi_from_bplus(T_VAR_Parameters *p) forecast : horizon x nvars matrix or null pointer horizon : positive integer - forecast horizon initial : initial value of predetermined variables. - shocks : array of length horizon of shocks or null pointer. If null - pointer, then the shocks are all zero. Each vector is of length + shocks : array of length horizon of shocks or null pointer. If null + pointer, then the shocks are all zero. Each vector is of length nvar. S : array of length horizon. S[t] is the state at time T+1+t. model : pointer to valid TStateModel structure. @@ -2408,7 +2408,7 @@ void Update_lambda_psi_from_bplus(T_VAR_Parameters *p) Computes forecast Returns: - The matrix forecast upon success and null upon failure. If forecast is + The matrix forecast upon success and null upon failure. If forecast is null, then it created. */ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *shocks, int *S, TStateModel *model) @@ -2455,7 +2455,7 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s } /* - For 1 <= k < h, y[k][i] is null if the ith coordinate of y(t0+1+k) is + For 1 <= k < h, y[k][i] is null if the ith coordinate of y(t0+1+k) is unrestricted and is its value otherwise. In general, t0 is the last index for which we have full information. It must be the case that t0 <= nobs. */ @@ -2570,7 +2570,7 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s /* } */ /* - + */ /* TVector* dw_state_space_mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model) */ /* { */ @@ -3580,7 +3580,7 @@ TMatrix MakeA0(TMatrix A0, int s, T_VAR_Parameters *p) /* Assumes: - A0 : Matrix array of length n_states or null pointer. A0[s] is either + A0 : Matrix array of length n_states or null pointer. A0[s] is either p->nvars x p->nvars matrix or null pointer */ TMatrix* MakeA0_All(TMatrix *A0, T_VAR_Parameters *p) @@ -3634,7 +3634,7 @@ TMatrix MakeAplus(TMatrix Aplus, int k, T_VAR_Parameters *p) /* Assumes: - Aplus : Matrix array of length n_states or null pointer. Aplus[s] is either + Aplus : Matrix array of length n_states or null pointer. Aplus[s] is either p->npre x p->nvars matrix or null pointer */ TMatrix* MakeAplus_All(TMatrix *Aplus, T_VAR_Parameters *p) @@ -3684,7 +3684,7 @@ TMatrix MakeZeta(TMatrix Zeta, int k, T_VAR_Parameters *p) /* Assumes: - Zeta : Matrix array of length n_states or null pointer. Zeta[s] is either + Zeta : Matrix array of length n_states or null pointer. Zeta[s] is either p->vars x p->nvars matrix or null pointer */ TMatrix* MakeZeta_All(TMatrix *Zeta, T_VAR_Parameters *p) diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index b590e5536..a332d33d9 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -20,7 +20,7 @@ model : point to valid TStateModel structure Results: - Computes and prints to the file f_out the requested percentiles for forecasts + Computes and prints to the file f_out the requested percentiles for forecasts of the observables. Returns: @@ -84,7 +84,7 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste { // Draw time T regime m=DrawDiscrete(init_prob); - + // Draw regimes from time T+1 through T+h inclusive for (j=0; j < h; j++) { @@ -138,7 +138,7 @@ ERROR_EXIT: model : point to valid TStateModel/T_MSStateSpace structure Results: - Computes and prints to the file f_out the requested percentiles for forecasts + Computes and prints to the file f_out the requested percentiles for forecasts of the observables. Returns: @@ -146,10 +146,10 @@ ERROR_EXIT: Notes: The regime at time T is drawn from the filtered probabilities at time t, and - is set to s there after. + is set to s there after. */ -int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, +int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, FILE *posterior_file, int s, int T, int h, TStateModel *model) { T_VAR_Parameters *p; @@ -238,54 +238,54 @@ ERROR_EXIT: /* - Attempt to set up model from command line. Command line options are the + Attempt to set up model from command line. Command line options are the following -ft If this argument exists, then the following is attempted: specification file name = est_final_.dat output file name = ir__regime_.dat - parameters file name = est_final_.dat + parameters file name = est_final_.dat header = "Posterior mode: " - + -fs - If this argument exists, then the specification file name is . + If this argument exists, then the specification file name is . The argument -fs takes precedence over -ft. -fp - If this argument exists, then the parameters file name is . The - argument -fp takes precedence over -ft. The default value is the filename + If this argument exists, then the parameters file name is . The + argument -fp takes precedence over -ft. The default value is the filename associated with the argument -fs. -ph
- If this argument exists, then the header for the parameters file is + If this argument exists, then the header for the parameters file is
. The default value is "Posterior mode: ". -horizon If this argument exists, then the horizon of the impulse responses is given by the passed integer. The default value is 12. - -error_bands + -error_bands Output error bands. (default = off - only median is computed) -percentiles n p_1 p_2 ... p_n - Percentiles to compute. The first parameter after percentiles must be the - number of percentiles and the following values are the actual percentiles. + Percentiles to compute. The first parameter after percentiles must be the + number of percentiles and the following values are the actual percentiles. default = 3 0.16 0.50 0.84 if error_bands flag is set = 1 0.50 otherwise - -parameter_uncertainty + -parameter_uncertainty Apply parameter uncertainty when computing error bands. - -shocks_per_parameter - Number of shocks and regime paths to draw for each parameter draw. The + -shocks_per_parameter + Number of shocks and regime paths to draw for each parameter draw. The default value is 1 if parameter_uncertainty is set and 10,000 otherwise. - - -thin - Thinning factor. Only 1/thin of the draws in posterior draws file are + + -thin + Thinning factor. Only 1/thin of the draws in posterior draws file are used. The default value is 1. - -regimes + -regimes Produces forecasts as if each regime were permanent. (default = off) -regime @@ -328,7 +328,7 @@ int main(int nargs, char **args) // parameter filename if (!parm) sprintf(parm=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); - } + } // horizon horizon=dw_ParseInteger_String(nargs,args,"horizon",12); @@ -344,7 +344,7 @@ int main(int nargs, char **args) " -horizon : horizon for the forecast (12)\n" ); exit(1); - } + } if (!parm) strcpy(parm=(char*)malloc(strlen(spec)+1),spec); @@ -364,7 +364,7 @@ int main(int nargs, char **args) free(head); free(parm); - //============================= Compute forecasts ============================= + //============================= Compute forecasts ============================= // Mean forecast /* if (dw_FindArgument_String(nargs,args,"mean") != -1) */ 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 6384e5948..be2c4b0c9 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -12,7 +12,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, +static TMatrix MakeHistogram(int low, int *h, PRECISION min, PRECISION max,int intervals, int sample_size, 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); @@ -53,7 +53,7 @@ TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int t if (!(h->high=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); for (i=rows-1; i >= 0; i--) - if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); h->Min=CreateMatrix(rows,cols); h->Max=CreateMatrix(rows,cols); @@ -476,8 +476,8 @@ static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int inte } /* - Adds a observation to the histogram. The type of the histogram must - be HISTOGRAM_VARIABLE. + Adds a observation to the histogram. The type of the histogram must + be HISTOGRAM_VARIABLE. */ static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) { @@ -485,7 +485,7 @@ static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISIO if ((x < *min) || (x > *max)) Resize(x,h,min,max,intervals); - if (*max > *min) + if (*max > *min) { i=(int)(intervals*(x - *min)/(*max - *min)); h[(i < intervals) ? i : intervals-1]++; @@ -516,13 +516,13 @@ static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECIS /* Returns the level such that the probability of observing an observation - less than or equal to level is percentile. If there is a point mass at + less than or equal to level is percentile. If there is a point mass at x, and P(y < x) <= percentile <= P(y <= x), then x is returned. Assumes Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also if - + non-negative. Also if + high = sample_size - (low + h[0] + ... + h[intervals - 1]), then high is non-negative. @@ -560,8 +560,8 @@ static PRECISION Percentile(PRECISION percentile, int low, int *h, PRECISION min Assumes Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also, if - + non-negative. Also, if + high = sample_size - (low + h[0] + ... + h[intervals - 1]), then high is non-negative. @@ -591,18 +591,18 @@ static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRE if (level >= max) return 1.0; for (count=low, i=0; i < intervals; count+=h[i++]) - if ((min+=inc) >= level) + if ((min+=inc) >= level) return ((PRECISION)count + (PRECISION)h[i]*(level - min + inc)/inc)/(PRECISION)sample_size; return 1.0; } /* Returns a histogram over the interval I=[min_out,max_out]. The matrix returned - has bins rows and 2 columns. If inc=(max_out - min_out)/bins, then the first - element of the ith row is + has bins rows and 2 columns. If inc=(max_out - min_out)/bins, then the first + element of the ith row is + + min + (i + 0.5)*inc, - min + (i + 0.5)*inc, - which is the mid-point of the ith interval. The second element is P(min + i*inc < x <= min + (i + 1)*inc)/inc, @@ -611,11 +611,11 @@ static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRE Assumes Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also if - + non-negative. Also if + high = sample_size - (low + h[0] + ... + h[intervals - 1]), - then high is non-negative. + then high is non-negative. If min < max, let inc=(max - min)/intervals and define @@ -742,7 +742,7 @@ static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECI } } } - + return MakeHistogram(low,h,min,max,intervals,sample_size,min_out,max_out,bins); } From 3bbfc05634b078e1fb30f198eb9a430d20244789 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 30 Aug 2010 15:01:26 +0200 Subject: [PATCH 04/12] SWZ: Change file ending type to unix --- .../DWCcode/histogram/dw_histogram.c | 1498 ++++++++--------- .../DWCcode/histogram/dw_histogram.h | 154 +- 2 files changed, 826 insertions(+), 826 deletions(-) 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 be2c4b0c9..79ef5bd4a 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -1,749 +1,749 @@ - -#include -#include -#include - -#include "dw_histogram.h" -#include "dw_error.h" - -static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); -static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); -static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals); - -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); -static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins); - - -/******************************************************************************* - The following set of routines create a matrix of histograms on the fly. -*******************************************************************************/ -/* - Assumes - rows > 0 - cols > 0 - intervals > 0 - type = HISTOGRAM_FIXED or HISTOGRAM_VARIABLE - - Results - Creates and returns a matrix histogram data structure. The size of the - matrix is m x n and the number of intervals is intrvls. - -*/ -TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type) -{ - int i, j; - TMatrixHistogram *h; - - if (!(h=(TMatrixHistogram *)malloc(sizeof(TMatrixHistogram)))) dw_Error(MEM_ERR); - - if (!(h->freq=(int***)malloc(rows*sizeof(int**)))) dw_Error(MEM_ERR); - for (i=rows-1; i >= 0; i--) - { - if (!(h->freq[i]=(int**)malloc(cols*sizeof(int*)))) dw_Error(MEM_ERR); - for (j=cols-1; j >= 0; j--) - if (!(h->freq[i][j]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); - } - - if (!(h->low=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); - for (i=rows-1; i >= 0; i--) - if (!(h->low[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); - - if (!(h->high=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); - for (i=rows-1; i >= 0; i--) - if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); - - h->Min=CreateMatrix(rows,cols); - h->Max=CreateMatrix(rows,cols); - - h->rows=rows; - h->cols=cols; - h->intervals=intervals; - h->sample_size=0; - h->type=type; - - return h; -} - -void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h) -{ - EquateMatrix(h->Min,Min); - EquateMatrix(h->Max,Max); - h->sample_size=0; -} - -void FreeMatrixHistogram(TMatrixHistogram *h) -{ - int i, j; - for (i=h->rows-1; i >= 0; i--) - { - for (j=h->cols-1; j >= 0; j--) free(h->freq[i][j]); - free(h->freq[i]); - } - free(h->freq); - for (i=h->rows-1; i >= 0; i--) free(h->low[i]); - free(h->low); - for (i=h->rows-1; i >= 0; i--) free(h->high[i]); - free(h->high); - FreeMatrix(h->Min); - FreeMatrix(h->Max); - free(h); -} - -void AddMatrixObservation(TMatrix X, TMatrixHistogram *h) -{ - int i, j, k; - - if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); - - if (h->sample_size <= 0) - { - 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; - } - 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); - } - - if (h->type == HISTOGRAM_FIXED) - for (i=h->rows-1; i >= 0; i--) - for (j=h->cols-1; j >= 0; j--) - AddObservationFixed(ElementM(X,i,j),h->low[i]+j,h->freq[i][j],h->high[i]+j,ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals); - else - for (i=h->rows-1; i >= 0; i--) - for (j=h->cols-1; j >= 0; j--) - AddObservationVariable(ElementM(X,i,j),h->freq[i][j],&ElementM(h->Min,i,j),&ElementM(h->Max,i,j),h->intervals); - - h->sample_size++; -} - -void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h) -{ - int i, j; - - if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); - - for (i=h->rows-1; i >= 0; i--) - for (j=h->cols-1; j >= 0; j--) - ElementM(X,i,j)=Percentile(percentile,h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); -} - -/* - Returns the probability that an observation is less than or equal to - level. - - Assumes - For 0 <= i < h->rows and 0 <= j < h->cols, let - - I[i][j][k]=(h->min[i][j] + k*inc[i][j], h->min[i][j] + (k+1)*inc[i][j]), - - where inc[i][j]=(h->max[i][j] - h->min[i][j])/h->samples_size. The - distribution is uniform on I[i][k][j] and - - P(h->min[i][j] + k*inc[i][j] < x[i][j] < h->min[i][j] + (k+1)*inc[i][j]) - = h->freq[i][j][k]/h->sample_size. - - Furthermore, - - P(x[i][j] < h->min[i][j]) = 0 and P(x[i][j] > h->min[i][j]) = 0. - - In addition, if h->type == FIXED, then - - P(x[i][j] = h->min[i][j]) = h->low[i][j]/h->sample_size - - and - - P(x[i][j] = h->min[i][j]) = h->high[i][j]/h->sample_size. -*/ -void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h) -{ - int i, j; - - if ((h->rows != RowM(P)) || (h->cols != ColM(P)) || - (h->rows != RowM(Level)) || (h->cols != ColM(Level))) - dw_Error(SIZE_ERR); - - for (i=h->rows-1; i >= 0; i--) - for (j=h->cols-1; j >= 0; j--) - ElementM(P,i,j)=Cumulative(ElementM(Level,i,j),h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); -} - -TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h) -{ - return MakeHistogramAuto(h->low[i][j],h->freq[i][j],h->high[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,bins); -} - -TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h) -{ - return MakeHistogram(h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,min,max,bins); -} - -/******************************************************************************* - The following set of routines create a vector of histograms on the fly. -*******************************************************************************/ -TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type) -{ - int i; - TVectorHistogram *h; - - if (!(h=(TVectorHistogram *)malloc(sizeof(TVectorHistogram)))) - dw_Error(MEM_ERR); - - if (!(h->freq=(int**)malloc(dim*sizeof(int*)))) dw_Error(MEM_ERR); - for (i=dim-1; i >= 0; i--) - if (!(h->freq[i]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); - - if (!(h->low=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); - if (!(h->high=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); - - h->Min=CreateVector(dim); - h->Max=CreateVector(dim); - - h->dim=dim; - h->intervals=intervals; - h->sample_size=0; - h->type=type; - - return h; -} - -void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h) -{ - EquateVector(h->Min,Min); - EquateVector(h->Max,Max); - h->sample_size=0; -} - -void FreeVectorHistogram(TVectorHistogram *h) -{ - int i; - for (i=h->dim-1; i >= 0; i--) free(h->freq[i]); - free(h->freq); - free(h->low); - free(h->high); - FreeVector(h->Min); - FreeVector(h->Max); - free(h); -} - -void AddVectorObservation(TVector x, TVectorHistogram *h) -{ - int i, k; - - if (h->dim != DimV(x)) dw_Error(SIZE_ERR); - - if (h->sample_size <= 0) - { - 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; - } - if (h->type == HISTOGRAM_VARIABLE) - for (i=h->dim-1; i >= 0; i--) - ElementV(h->Min,i)=ElementV(h->Max,i)=ElementV(x,i); - } - - if (h->type == HISTOGRAM_FIXED) - for (i=h->dim-1; i >= 0; i--) - AddObservationFixed(ElementV(x,i),h->low+i,h->freq[i],h->high+i,ElementV(h->Min,i),ElementV(h->Max,i),h->intervals); - else - for (i=h->dim-1; i >= 0; i--) - AddObservationVariable(ElementV(x,i),h->freq[i],&ElementV(h->Min,i),&ElementV(h->Max,i),h->intervals); - - h->sample_size++; -} - -void VectorPercentile(TVector x, PRECISION percentile, TVectorHistogram *h) -{ - int i; - - if (h->dim != DimV(x)) dw_Error(SIZE_ERR); - - for (i=h->dim-1; i >= 0; i--) - ElementV(x,i)=Percentile(percentile,h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); -} - - - -/* - Returns the probability that an observation is less than or equal to - level. - - Assumes - For 0 <= i < h->dim, let - - I[i][k]=(h->min[i] + k*inc[i], h->min[i] + (k+1)*inc[i]), - - where inc[i]=(h->max[i] - h->min[i])/h->samples_size. The distribution - is uniform on I[i][k] and - - P(h->min[i] + k*inc[i] < x[i] < h->min[i] + (k+1)*inc[i]) - = h->freq[i][k]/h->sample_size. - - Furthermore, - - P(x[i] < h->min[i]) = 0 and P(x[i] > h->min[i]) = 0. - - In addition, if h->type == FIXED, then - - P(x[i] = h->min[i]) = h->low[i]/h->sample_size - - and - - P(x[i] = h->min[i]) = h->high[i]/h->sample_size. -*/ -void VectorCumulative(TVector p, TVector level, TVectorHistogram *h) -{ - int i; - - if (h->dim != DimV(p) || (h->dim != DimV(level))) - dw_Error(SIZE_ERR); - - for (i=h->dim-1; i >= 0; i--) - ElementV(p,i)=Cumulative(ElementV(level,i),h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); -} - -TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h) -{ - return MakeHistogramAuto(h->low[i],h->freq[i],h->high[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,bins); -} - -TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h) -{ - return MakeHistogram(h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,min,max,bins); -} -/******************************************************************************* - The following set of routines create a scalar histogram on the fly. -*******************************************************************************/ -/* - Assumes - - Results - Creates and returns a scalar histogram data structure. -*/ -TScalarHistogram *CreateScalarHistogram(int intervals, int type) -{ - TScalarHistogram *h; - - if (!(h=(TScalarHistogram *)malloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR); - - if (!(h->freq=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); - - h->intervals=intervals; - h->sample_size=0; - h->type=type; - - return h; -} - -void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h) -{ - h->Min=Min; - h->Max=Max; - h->sample_size=0; -} - -void FreeScalarHistogram(TScalarHistogram *h) -{ - free(h->freq); - free(h); -} - -void AddScalarObservation(PRECISION x, TScalarHistogram *h) -{ - int k; - - if (h->sample_size <= 0) - { - h->low=h->high=0; - for (k=h->intervals-1; k >= 0; k--) h->freq[k]=0; - if (h->type == HISTOGRAM_VARIABLE) h->Min=h->Max=x; - } - - if (h->type == HISTOGRAM_FIXED) - AddObservationFixed(x,&(h->low),h->freq,&(h->high),h->Min,h->Max,h->intervals); - else - AddObservationVariable(x,h->freq,&(h->Min),&(h->Max),h->intervals); - - h->sample_size++; -} - -PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h) -{ - return Percentile(percentile,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); -} - -/* - Returns the probability that an observation is less than or equal to - level. - - Assumes - Let - - I[k]=(h->min + k*inc, h->min + (k+1)*inc), - - where inc=(h->max - h->min)/h->samples_size. The distribution - is uniform on I[k] and - - P(h->min + k*inc < x < h->min + (k+1)*inc) = h->freq[k]/h->sample_size. - - Furthermore, - - P(x < h->min) = 0 and P(x > h->min) = 0. - - In addition, if h->type == FIXED, then - - P(x = h->min) = h->low/h->sample_size - - and - - P(x = h->min) = h->high/h->sample_size. -*/ -PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h) -{ - return Cumulative(level,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); -} - -TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h) -{ - return MakeHistogramAuto(h->low,h->freq,h->Min,h->high,h->Max,h->intervals,h->sample_size,bins); -} - -TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h) -{ - return MakeHistogram(h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size,min,max,bins); -} - -/*******************************************************************************/ -/***************************** Low Level Routines ******************************/ -/*******************************************************************************/ -/* - Resizes the histogram. After resizing, it is guaranteed that *min <= x <= *max. - The type of the histogram must be HISTOGRAM_VARIABLE. -*/ -static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) -{ - int i, j, k, m; - if (x > *max) - if (x - *min >= (PRECISION)intervals*(*max - *min)) - { - for (i=1; i < intervals; i++) - { - h[0]+=h[i]; - h[i]=0; - } - *max=x; - } - else - { - m=(int)ceil((x - *min)/(*max - *min)); - for (i=j=0; i < intervals; j++) - for(h[j]=h[i++], k=1; (k < m) && (i < intervals); k++) - h[j]+=h[i++]; - for ( ; j < intervals; j++) h[j]=0; - *max=*min + m*(*max - *min); - if (x > *max) *max=x; - } - else - if (x < *min) - if (*max - x >= (PRECISION)intervals*(*max - *min)) - { - for (j=intervals-1, i=intervals-2; i >= 0; i--) - { - h[j]+=h[i]; - h[i]=0; - } - *min=x; - } - else - { - m=(int)ceil((*max - x)/(*max - *min)); - for (i=j=intervals-1; i >= 0; j--) - for(h[j]=h[i--], k=1; (k < m) && (i >= 0); k++) - h[j]+=h[i--]; - for ( ; j >= 0; j--) h[j]=0; - *min=*max - m*(*max - *min); - if (x < *min) *min=x; - } -} - -/* - Adds a observation to the histogram. The type of the histogram must - be HISTOGRAM_VARIABLE. -*/ -static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) -{ - int i; - - if ((x < *min) || (x > *max)) Resize(x,h,min,max,intervals); - - if (*max > *min) - { - i=(int)(intervals*(x - *min)/(*max - *min)); - h[(i < intervals) ? i : intervals-1]++; - } - else - h[0]++; -} - -/* - Adds a observation to the histogram. The type of the histogram must - be HISTOGRAM_FIXED. -*/ -static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals) -{ - PRECISION y=floor(intervals*(x - min)/(max - min)); - if (y < 0) - (*low)++; - else - if (y < intervals) - h[(int)y]++; - else - (*high)++; -} - -/******************************************************************************/ -/******************************************************************************/ -/******************************************************************************/ - -/* - Returns the level such that the probability of observing an observation - less than or equal to level is percentile. If there is a point mass at - x, and P(y < x) <= percentile <= P(y <= x), then x is returned. - - Assumes - Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also if - - high = sample_size - (low + h[0] + ... + h[intervals - 1]), - - then high is non-negative. - - If min < max, let inc=(max - min)/intervals and define - - I[k]=(min + k*inc, min + (k+1)*inc), - - The distribution is uniform on I[k] and - - P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. - - Furthermore, there are point masses at min and max with probability - - P(x = min) = low/sample_size - and - P(x = max) = high/sample_size. - - If min = max, then there is a single point mass at this point. -*/ -static PRECISION Percentile(PRECISION percentile, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) -{ - int i; - percentile=percentile*sample_size - low; - if (percentile <= 0) return min; - for (i=0; i < intervals; i++) - if (h[i] && (percentile-=h[i]) <= 0) - return min + ((PRECISION)(i+1) + percentile/(PRECISION)h[i])*(max - min)/(PRECISION)intervals; - return max; -} - -/* - Returns the probability that an observation is less than or equal to - level. - - Assumes - Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also, if - - high = sample_size - (low + h[0] + ... + h[intervals - 1]), - - then high is non-negative. - - If min < max, let inc=(max - min)/intervals and define - - I[k]=(min + k*inc, min + (k+1)*inc), - - The distribution is uniform on I[k] and - - P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. - - Furthermore, there are point masses at min and max with probability - - P(x = min) = low/sample_size - and - P(x = max) = high/sample_size. - - If min = max, then there is a single point mass at this point -*/ -static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) -{ - PRECISION inc=(max-min)/(PRECISION)intervals; - int i, count; - - if (level < min) return 0.0; - if (level >= max) return 1.0; - - for (count=low, i=0; i < intervals; count+=h[i++]) - if ((min+=inc) >= level) - return ((PRECISION)count + (PRECISION)h[i]*(level - min + inc)/inc)/(PRECISION)sample_size; - return 1.0; -} - -/* - Returns a histogram over the interval I=[min_out,max_out]. The matrix returned - has bins rows and 2 columns. If inc=(max_out - min_out)/bins, then the first - element of the ith row is - - min + (i + 0.5)*inc, - - which is the mid-point of the ith interval. The second element is - - P(min + i*inc < x <= min + (i + 1)*inc)/inc, - - which is the average density over the ith interval. - - Assumes - Both intervals and sample_size are poitive and low and h[i] are - non-negative. Also if - - high = sample_size - (low + h[0] + ... + h[intervals - 1]), - - then high is non-negative. - - If min < max, let inc=(max - min)/intervals and define - - I[k]=(min + k*inc, min + (k+1)*inc), - - The distribution is uniform on I[k] and - - P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. - - Furthermore, there are point masses at min and max with probability - - P(x = min) = low/sample_size - and - P(x = max) = high/sample_size. - - If min = max, then there is a single point mass at this point. -*/ -static TMatrix MakeHistogram(int low, int *h, PRECISION min, PRECISION max, - int intervals, int sample_size, PRECISION min_out, PRECISION max_out, int bins) -{ - int i; - PRECISION inc, x, cdf_lower, cdf_upper; - TMatrix X; - - inc=(max_out-min_out)/(PRECISION)bins; - - if (inc > 0) - { - X=CreateMatrix(bins,2); - x=min_out+inc; - - cdf_lower=Cumulative(min_out,low,h,min,max,intervals,sample_size); - - for (i=0; i < bins; i++) - { - cdf_upper=Cumulative(x,low,h,min,max,intervals,sample_size); - - ElementM(X,i,0)=x - 0.5*inc; - ElementM(X,i,1)=(cdf_upper-cdf_lower)/inc; - - cdf_lower=cdf_upper; - x+=inc; - } - } - else - return (TMatrix)NULL; - - return X; -} - -/* - Automatically chooses lenth of interval over which to produce histogram and - then calls MakeHistogram(). -*/ -static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins) -{ - PRECISION inc=(max-min)/intervals, max_out, min_out; - int lo, hi; - - if ((low == sample_size) || (inc <= 0)) - { - min_out=min-1.0; - max_out=min+1.0; - } - else - { - if (low > 0) - lo=-1; - else - for (lo=0; (lo < intervals) && !h[lo]; lo++); - - if (lo == intervals) - { - 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 (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; - } - } - } - - return MakeHistogram(low,h,min,max,intervals,sample_size,min_out,max_out,bins); -} - - + +#include +#include +#include + +#include "dw_histogram.h" +#include "dw_error.h" + +static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); +static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); +static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals); + +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); +static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins); + + +/******************************************************************************* + The following set of routines create a matrix of histograms on the fly. +*******************************************************************************/ +/* + Assumes + rows > 0 + cols > 0 + intervals > 0 + type = HISTOGRAM_FIXED or HISTOGRAM_VARIABLE + + Results + Creates and returns a matrix histogram data structure. The size of the + matrix is m x n and the number of intervals is intrvls. + +*/ +TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type) +{ + int i, j; + TMatrixHistogram *h; + + if (!(h=(TMatrixHistogram *)malloc(sizeof(TMatrixHistogram)))) dw_Error(MEM_ERR); + + if (!(h->freq=(int***)malloc(rows*sizeof(int**)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + { + if (!(h->freq[i]=(int**)malloc(cols*sizeof(int*)))) dw_Error(MEM_ERR); + for (j=cols-1; j >= 0; j--) + if (!(h->freq[i][j]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + } + + if (!(h->low=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + if (!(h->low[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + + if (!(h->high=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=rows-1; i >= 0; i--) + if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + + h->Min=CreateMatrix(rows,cols); + h->Max=CreateMatrix(rows,cols); + + h->rows=rows; + h->cols=cols; + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h) +{ + EquateMatrix(h->Min,Min); + EquateMatrix(h->Max,Max); + h->sample_size=0; +} + +void FreeMatrixHistogram(TMatrixHistogram *h) +{ + int i, j; + for (i=h->rows-1; i >= 0; i--) + { + for (j=h->cols-1; j >= 0; j--) free(h->freq[i][j]); + free(h->freq[i]); + } + free(h->freq); + for (i=h->rows-1; i >= 0; i--) free(h->low[i]); + free(h->low); + for (i=h->rows-1; i >= 0; i--) free(h->high[i]); + free(h->high); + FreeMatrix(h->Min); + FreeMatrix(h->Max); + free(h); +} + +void AddMatrixObservation(TMatrix X, TMatrixHistogram *h) +{ + int i, j, k; + + if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); + + if (h->sample_size <= 0) + { + 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; + } + 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); + } + + if (h->type == HISTOGRAM_FIXED) + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + AddObservationFixed(ElementM(X,i,j),h->low[i]+j,h->freq[i][j],h->high[i]+j,ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals); + else + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + AddObservationVariable(ElementM(X,i,j),h->freq[i][j],&ElementM(h->Min,i,j),&ElementM(h->Max,i,j),h->intervals); + + h->sample_size++; +} + +void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h) +{ + int i, j; + + if ((h->rows != RowM(X)) || (h->cols != ColM(X))) dw_Error(SIZE_ERR); + + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + ElementM(X,i,j)=Percentile(percentile,h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + For 0 <= i < h->rows and 0 <= j < h->cols, let + + I[i][j][k]=(h->min[i][j] + k*inc[i][j], h->min[i][j] + (k+1)*inc[i][j]), + + where inc[i][j]=(h->max[i][j] - h->min[i][j])/h->samples_size. The + distribution is uniform on I[i][k][j] and + + P(h->min[i][j] + k*inc[i][j] < x[i][j] < h->min[i][j] + (k+1)*inc[i][j]) + = h->freq[i][j][k]/h->sample_size. + + Furthermore, + + P(x[i][j] < h->min[i][j]) = 0 and P(x[i][j] > h->min[i][j]) = 0. + + In addition, if h->type == FIXED, then + + P(x[i][j] = h->min[i][j]) = h->low[i][j]/h->sample_size + + and + + P(x[i][j] = h->min[i][j]) = h->high[i][j]/h->sample_size. +*/ +void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h) +{ + int i, j; + + if ((h->rows != RowM(P)) || (h->cols != ColM(P)) || + (h->rows != RowM(Level)) || (h->cols != ColM(Level))) + dw_Error(SIZE_ERR); + + for (i=h->rows-1; i >= 0; i--) + for (j=h->cols-1; j >= 0; j--) + ElementM(P,i,j)=Cumulative(ElementM(Level,i,j),h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size); +} + +TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h) +{ + return MakeHistogramAuto(h->low[i][j],h->freq[i][j],h->high[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,bins); +} + +TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h) +{ + return MakeHistogram(h->low[i][j],h->freq[i][j],ElementM(h->Min,i,j),ElementM(h->Max,i,j),h->intervals,h->sample_size,min,max,bins); +} + +/******************************************************************************* + The following set of routines create a vector of histograms on the fly. +*******************************************************************************/ +TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type) +{ + int i; + TVectorHistogram *h; + + if (!(h=(TVectorHistogram *)malloc(sizeof(TVectorHistogram)))) + dw_Error(MEM_ERR); + + if (!(h->freq=(int**)malloc(dim*sizeof(int*)))) dw_Error(MEM_ERR); + for (i=dim-1; i >= 0; i--) + if (!(h->freq[i]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + + if (!(h->low=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->high=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + + h->Min=CreateVector(dim); + h->Max=CreateVector(dim); + + h->dim=dim; + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h) +{ + EquateVector(h->Min,Min); + EquateVector(h->Max,Max); + h->sample_size=0; +} + +void FreeVectorHistogram(TVectorHistogram *h) +{ + int i; + for (i=h->dim-1; i >= 0; i--) free(h->freq[i]); + free(h->freq); + free(h->low); + free(h->high); + FreeVector(h->Min); + FreeVector(h->Max); + free(h); +} + +void AddVectorObservation(TVector x, TVectorHistogram *h) +{ + int i, k; + + if (h->dim != DimV(x)) dw_Error(SIZE_ERR); + + if (h->sample_size <= 0) + { + 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; + } + if (h->type == HISTOGRAM_VARIABLE) + for (i=h->dim-1; i >= 0; i--) + ElementV(h->Min,i)=ElementV(h->Max,i)=ElementV(x,i); + } + + if (h->type == HISTOGRAM_FIXED) + for (i=h->dim-1; i >= 0; i--) + AddObservationFixed(ElementV(x,i),h->low+i,h->freq[i],h->high+i,ElementV(h->Min,i),ElementV(h->Max,i),h->intervals); + else + for (i=h->dim-1; i >= 0; i--) + AddObservationVariable(ElementV(x,i),h->freq[i],&ElementV(h->Min,i),&ElementV(h->Max,i),h->intervals); + + h->sample_size++; +} + +void VectorPercentile(TVector x, PRECISION percentile, TVectorHistogram *h) +{ + int i; + + if (h->dim != DimV(x)) dw_Error(SIZE_ERR); + + for (i=h->dim-1; i >= 0; i--) + ElementV(x,i)=Percentile(percentile,h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); +} + + + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + For 0 <= i < h->dim, let + + I[i][k]=(h->min[i] + k*inc[i], h->min[i] + (k+1)*inc[i]), + + where inc[i]=(h->max[i] - h->min[i])/h->samples_size. The distribution + is uniform on I[i][k] and + + P(h->min[i] + k*inc[i] < x[i] < h->min[i] + (k+1)*inc[i]) + = h->freq[i][k]/h->sample_size. + + Furthermore, + + P(x[i] < h->min[i]) = 0 and P(x[i] > h->min[i]) = 0. + + In addition, if h->type == FIXED, then + + P(x[i] = h->min[i]) = h->low[i]/h->sample_size + + and + + P(x[i] = h->min[i]) = h->high[i]/h->sample_size. +*/ +void VectorCumulative(TVector p, TVector level, TVectorHistogram *h) +{ + int i; + + if (h->dim != DimV(p) || (h->dim != DimV(level))) + dw_Error(SIZE_ERR); + + for (i=h->dim-1; i >= 0; i--) + ElementV(p,i)=Cumulative(ElementV(level,i),h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size); +} + +TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h) +{ + return MakeHistogramAuto(h->low[i],h->freq[i],h->high[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,bins); +} + +TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h) +{ + return MakeHistogram(h->low[i],h->freq[i],ElementV(h->Min,i),ElementV(h->Max,i),h->intervals,h->sample_size,min,max,bins); +} +/******************************************************************************* + The following set of routines create a scalar histogram on the fly. +*******************************************************************************/ +/* + Assumes + + Results + Creates and returns a scalar histogram data structure. +*/ +TScalarHistogram *CreateScalarHistogram(int intervals, int type) +{ + TScalarHistogram *h; + + if (!(h=(TScalarHistogram *)malloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR); + + if (!(h->freq=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + + h->intervals=intervals; + h->sample_size=0; + h->type=type; + + return h; +} + +void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h) +{ + h->Min=Min; + h->Max=Max; + h->sample_size=0; +} + +void FreeScalarHistogram(TScalarHistogram *h) +{ + free(h->freq); + free(h); +} + +void AddScalarObservation(PRECISION x, TScalarHistogram *h) +{ + int k; + + if (h->sample_size <= 0) + { + h->low=h->high=0; + for (k=h->intervals-1; k >= 0; k--) h->freq[k]=0; + if (h->type == HISTOGRAM_VARIABLE) h->Min=h->Max=x; + } + + if (h->type == HISTOGRAM_FIXED) + AddObservationFixed(x,&(h->low),h->freq,&(h->high),h->Min,h->Max,h->intervals); + else + AddObservationVariable(x,h->freq,&(h->Min),&(h->Max),h->intervals); + + h->sample_size++; +} + +PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h) +{ + return Percentile(percentile,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + Let + + I[k]=(h->min + k*inc, h->min + (k+1)*inc), + + where inc=(h->max - h->min)/h->samples_size. The distribution + is uniform on I[k] and + + P(h->min + k*inc < x < h->min + (k+1)*inc) = h->freq[k]/h->sample_size. + + Furthermore, + + P(x < h->min) = 0 and P(x > h->min) = 0. + + In addition, if h->type == FIXED, then + + P(x = h->min) = h->low/h->sample_size + + and + + P(x = h->min) = h->high/h->sample_size. +*/ +PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h) +{ + return Cumulative(level,h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size); +} + +TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h) +{ + return MakeHistogramAuto(h->low,h->freq,h->Min,h->high,h->Max,h->intervals,h->sample_size,bins); +} + +TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h) +{ + return MakeHistogram(h->low,h->freq,h->Min,h->Max,h->intervals,h->sample_size,min,max,bins); +} + +/*******************************************************************************/ +/***************************** Low Level Routines ******************************/ +/*******************************************************************************/ +/* + Resizes the histogram. After resizing, it is guaranteed that *min <= x <= *max. + The type of the histogram must be HISTOGRAM_VARIABLE. +*/ +static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) +{ + int i, j, k, m; + if (x > *max) + if (x - *min >= (PRECISION)intervals*(*max - *min)) + { + for (i=1; i < intervals; i++) + { + h[0]+=h[i]; + h[i]=0; + } + *max=x; + } + else + { + m=(int)ceil((x - *min)/(*max - *min)); + for (i=j=0; i < intervals; j++) + for(h[j]=h[i++], k=1; (k < m) && (i < intervals); k++) + h[j]+=h[i++]; + for ( ; j < intervals; j++) h[j]=0; + *max=*min + m*(*max - *min); + if (x > *max) *max=x; + } + else + if (x < *min) + if (*max - x >= (PRECISION)intervals*(*max - *min)) + { + for (j=intervals-1, i=intervals-2; i >= 0; i--) + { + h[j]+=h[i]; + h[i]=0; + } + *min=x; + } + else + { + m=(int)ceil((*max - x)/(*max - *min)); + for (i=j=intervals-1; i >= 0; j--) + for(h[j]=h[i--], k=1; (k < m) && (i >= 0); k++) + h[j]+=h[i--]; + for ( ; j >= 0; j--) h[j]=0; + *min=*max - m*(*max - *min); + if (x < *min) *min=x; + } +} + +/* + Adds a observation to the histogram. The type of the histogram must + be HISTOGRAM_VARIABLE. +*/ +static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals) +{ + int i; + + if ((x < *min) || (x > *max)) Resize(x,h,min,max,intervals); + + if (*max > *min) + { + i=(int)(intervals*(x - *min)/(*max - *min)); + h[(i < intervals) ? i : intervals-1]++; + } + else + h[0]++; +} + +/* + Adds a observation to the histogram. The type of the histogram must + be HISTOGRAM_FIXED. +*/ +static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals) +{ + PRECISION y=floor(intervals*(x - min)/(max - min)); + if (y < 0) + (*low)++; + else + if (y < intervals) + h[(int)y]++; + else + (*high)++; +} + +/******************************************************************************/ +/******************************************************************************/ +/******************************************************************************/ + +/* + Returns the level such that the probability of observing an observation + less than or equal to level is percentile. If there is a point mass at + x, and P(y < x) <= percentile <= P(y <= x), then x is returned. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point. +*/ +static PRECISION Percentile(PRECISION percentile, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) +{ + int i; + percentile=percentile*sample_size - low; + if (percentile <= 0) return min; + for (i=0; i < intervals; i++) + if (h[i] && (percentile-=h[i]) <= 0) + return min + ((PRECISION)(i+1) + percentile/(PRECISION)h[i])*(max - min)/(PRECISION)intervals; + return max; +} + +/* + Returns the probability that an observation is less than or equal to + level. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also, if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point +*/ +static PRECISION Cumulative(PRECISION level, int low, int *h, PRECISION min, PRECISION max, int intervals, int sample_size) +{ + PRECISION inc=(max-min)/(PRECISION)intervals; + int i, count; + + if (level < min) return 0.0; + if (level >= max) return 1.0; + + for (count=low, i=0; i < intervals; count+=h[i++]) + if ((min+=inc) >= level) + return ((PRECISION)count + (PRECISION)h[i]*(level - min + inc)/inc)/(PRECISION)sample_size; + return 1.0; +} + +/* + Returns a histogram over the interval I=[min_out,max_out]. The matrix returned + has bins rows and 2 columns. If inc=(max_out - min_out)/bins, then the first + element of the ith row is + + min + (i + 0.5)*inc, + + which is the mid-point of the ith interval. The second element is + + P(min + i*inc < x <= min + (i + 1)*inc)/inc, + + which is the average density over the ith interval. + + Assumes + Both intervals and sample_size are poitive and low and h[i] are + non-negative. Also if + + high = sample_size - (low + h[0] + ... + h[intervals - 1]), + + then high is non-negative. + + If min < max, let inc=(max - min)/intervals and define + + I[k]=(min + k*inc, min + (k+1)*inc), + + The distribution is uniform on I[k] and + + P(min + k*inc < x < min + (k+1)*inc) = h[k]/sample_size. + + Furthermore, there are point masses at min and max with probability + + P(x = min) = low/sample_size + and + P(x = max) = high/sample_size. + + If min = max, then there is a single point mass at this point. +*/ +static TMatrix MakeHistogram(int low, int *h, PRECISION min, PRECISION max, + int intervals, int sample_size, PRECISION min_out, PRECISION max_out, int bins) +{ + int i; + PRECISION inc, x, cdf_lower, cdf_upper; + TMatrix X; + + inc=(max_out-min_out)/(PRECISION)bins; + + if (inc > 0) + { + X=CreateMatrix(bins,2); + x=min_out+inc; + + cdf_lower=Cumulative(min_out,low,h,min,max,intervals,sample_size); + + for (i=0; i < bins; i++) + { + cdf_upper=Cumulative(x,low,h,min,max,intervals,sample_size); + + ElementM(X,i,0)=x - 0.5*inc; + ElementM(X,i,1)=(cdf_upper-cdf_lower)/inc; + + cdf_lower=cdf_upper; + x+=inc; + } + } + else + return (TMatrix)NULL; + + return X; +} + +/* + Automatically chooses lenth of interval over which to produce histogram and + then calls MakeHistogram(). +*/ +static TMatrix MakeHistogramAuto(int low, int *h, int high, PRECISION min, PRECISION max, int intervals, int sample_size, int bins) +{ + PRECISION inc=(max-min)/intervals, max_out, min_out; + int lo, hi; + + if ((low == sample_size) || (inc <= 0)) + { + min_out=min-1.0; + max_out=min+1.0; + } + else + { + if (low > 0) + lo=-1; + else + for (lo=0; (lo < intervals) && !h[lo]; lo++); + + if (lo == intervals) + { + 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 (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; + } + } + } + + return MakeHistogram(low,h,min,max,intervals,sample_size,min_out,max_out,bins); +} + + diff --git a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h index 801a8859d..64e24a2e3 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h @@ -1,77 +1,77 @@ -#ifndef __HISTOGRAMS__ -#define __HISTOGRAMS__ - -#include "matrix.h" - -#define HISTOGRAM_FIXED 1 -#define HISTOGRAM_VARIABLE 2 - -/* Matrix histograms */ -typedef struct -{ - TMatrix Min; - TMatrix Max; - int **low; - int **high; - int ***freq; - int rows; - int cols; - int intervals; - int sample_size; - int type; -} TMatrixHistogram; - -/* Vector histograms */ -typedef struct -{ - TVector Min; - TVector Max; - int *low; - int *high; - int **freq; - int dim; - int intervals; - int sample_size; - int type; -} TVectorHistogram; - -/* Scalar histograms */ -typedef struct -{ - PRECISION Min; - PRECISION Max; - int low; - int high; - int *freq; - int intervals; - int sample_size; - int type; -} TScalarHistogram; - -TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type); -void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h); -void FreeMatrixHistogram(TMatrixHistogram *h); -void AddMatrixObservation(TMatrix X, TMatrixHistogram *h); -void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h); -void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h); -TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h); -TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h); - -TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type); -void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h); -void FreeVectorHistogram(TVectorHistogram *h); -void AddVectorObservation(TVector X, TVectorHistogram *h); -void VectorPercentile(TVector X, PRECISION percentile, TVectorHistogram *h); -void VectorCumulative(TVector p, TVector level, TVectorHistogram *h); -TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h); -TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h); - -TScalarHistogram *CreateScalarHistogram(int intervals, int type); -void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h); -void FreeScalarHistogram(TScalarHistogram *h); -void AddScalarObservation(PRECISION x, TScalarHistogram *h); -PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h); -PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h); -TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h); -TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h); -#endif +#ifndef __HISTOGRAMS__ +#define __HISTOGRAMS__ + +#include "matrix.h" + +#define HISTOGRAM_FIXED 1 +#define HISTOGRAM_VARIABLE 2 + +/* Matrix histograms */ +typedef struct +{ + TMatrix Min; + TMatrix Max; + int **low; + int **high; + int ***freq; + int rows; + int cols; + int intervals; + int sample_size; + int type; +} TMatrixHistogram; + +/* Vector histograms */ +typedef struct +{ + TVector Min; + TVector Max; + int *low; + int *high; + int **freq; + int dim; + int intervals; + int sample_size; + int type; +} TVectorHistogram; + +/* Scalar histograms */ +typedef struct +{ + PRECISION Min; + PRECISION Max; + int low; + int high; + int *freq; + int intervals; + int sample_size; + int type; +} TScalarHistogram; + +TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int type); +void SetMaxMinMatrixHistogram(TMatrix Min, TMatrix Max, TMatrixHistogram *h); +void FreeMatrixHistogram(TMatrixHistogram *h); +void AddMatrixObservation(TMatrix X, TMatrixHistogram *h); +void MatrixPercentile(TMatrix X, PRECISION percentile, TMatrixHistogram *h); +void MatrixCumulative(TMatrix P, TMatrix Level, TMatrixHistogram *h); +TMatrix PlotMatrixHistogramAuto(int i, int j, int bins, TMatrixHistogram *h); +TMatrix PlotMatrixHistogram(int i, int j, PRECISION min, PRECISION max, int bins, TMatrixHistogram *h); + +TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type); +void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h); +void FreeVectorHistogram(TVectorHistogram *h); +void AddVectorObservation(TVector X, TVectorHistogram *h); +void VectorPercentile(TVector X, PRECISION percentile, TVectorHistogram *h); +void VectorCumulative(TVector p, TVector level, TVectorHistogram *h); +TMatrix PlotVectorHistogramAuto(int i, int bins, TVectorHistogram *h); +TMatrix PlotVectorHistogram(int i, PRECISION min, PRECISION max, int bins, TVectorHistogram *h); + +TScalarHistogram *CreateScalarHistogram(int intervals, int type); +void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h); +void FreeScalarHistogram(TScalarHistogram *h); +void AddScalarObservation(PRECISION x, TScalarHistogram *h); +PRECISION ScalarPercentile(PRECISION percentile, TScalarHistogram *h); +PRECISION ScalarCumulative(PRECISION level, TScalarHistogram *h); +TMatrix PlotScalarHistogramAuto(int bins, TScalarHistogram *h); +TMatrix PlotScalarHistogram(PRECISION min, PRECISION max, int bins, TScalarHistogram *h); +#endif From 12993a810883740f19e7fa7f37502ca472c7dee8 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 17:58:56 +0200 Subject: [PATCH 05/12] SWZ: replace malloc with swzMalloc --- matlab/swz/c-code/sbvar/var/forecast.c | 30 +++++++++-------- .../DWCcode/histogram/dw_histogram.c | 32 ++++++++++--------- 2 files changed, 33 insertions(+), 29 deletions(-) diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index a332d33d9..d98baebd0 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -9,6 +9,8 @@ #include #include +#include "modify_for_mex.h" + /* Assumes f_out : valid FILE pointer @@ -43,7 +45,7 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste if (T > p->nobs) return 0; // allocate memory - S=(int*)malloc(h*sizeof(int)); + S=(int*)swzMalloc(h*sizeof(int)); forecast=CreateMatrix(h,p->nvars); histogram=CreateMatrixHistogram(h,p->nvars,100,HISTOGRAM_VARIABLE); initial=CreateVector(p->npre); @@ -166,7 +168,7 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, if (T > p->nobs) return 0; // allocate memory - S=(int*)malloc(h*sizeof(int)); + S=(int*)swzMalloc(h*sizeof(int)); for (i=0; i < h; i++) S[i]=s; forecast=CreateMatrix(h,p->nvars); histogram=CreateMatrixHistogram(h,p->nvars,100,HISTOGRAM_VARIABLE); @@ -306,15 +308,15 @@ int main(int nargs, char **args) // specification filename if (buffer=dw_ParseString_String(nargs,args,"fs",(char*)NULL)) - strcpy(spec=(char*)malloc(strlen(buffer)+1),buffer); + strcpy(spec=(char*)swzMalloc(strlen(buffer)+1),buffer); // parameter filename if (buffer=dw_ParseString_String(nargs,args,"fp",(char*)NULL)) - strcpy(parm=(char*)malloc(strlen(buffer)+1),buffer); + strcpy(parm=(char*)swzMalloc(strlen(buffer)+1),buffer); // header if (buffer=dw_ParseString_String(nargs,args,"ph",(char*)NULL)) - strcpy(head=(char*)malloc(strlen(buffer)+1),buffer); + strcpy(head=(char*)swzMalloc(strlen(buffer)+1),buffer); // file tag if (tag=dw_ParseString_String(nargs,args,"ft",(char*)NULL)) @@ -323,11 +325,11 @@ int main(int nargs, char **args) // specification filename if (!spec) - sprintf(spec=(char*)malloc(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*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + sprintf(parm=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); } // horizon @@ -347,12 +349,12 @@ int main(int nargs, char **args) } if (!parm) - strcpy(parm=(char*)malloc(strlen(spec)+1),spec); + strcpy(parm=(char*)swzMalloc(strlen(spec)+1),spec); if (!head) { buffer="Posterior mode: "; - strcpy(head=(char*)malloc(strlen(buffer)+1),buffer); + strcpy(head=(char*)swzMalloc(strlen(buffer)+1),buffer); } model=Read_VAR_Specification((FILE*)NULL,spec); @@ -370,7 +372,7 @@ int main(int nargs, char **args) /* if (dw_FindArgument_String(nargs,args,"mean") != -1) */ /* { */ /* fmt="forecasts_mean_%s.prn"; */ - /* sprintf(out_filename=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); */ + /* sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); */ /* f_out=fopen(out_filename,"wt"); */ /* free(out_filename); */ /* printf("Constructing mean forecast\n"); */ @@ -386,7 +388,7 @@ int main(int nargs, char **args) { // Open posterior draws file fmt="draws_%s.dat"; - sprintf(post=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + 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); @@ -448,7 +450,7 @@ int main(int nargs, char **args) { rewind(posterior_file); fmt="forecasts_percentiles_regime_%d_%s.prn"; - sprintf(out_filename=(char*)malloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); f_out=fopen(out_filename,"wt"); free(out_filename); printf("Constructing percentiles for forecasts - regime %d\n",s); @@ -459,7 +461,7 @@ int main(int nargs, char **args) if (((s=dw_ParseInteger_String(nargs,args,"regime",-1)) >= 0) && (s < p->nstates)) { fmt="forecasts_percentiles_regime_%d_%s.prn"; - sprintf(out_filename=(char*)malloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 3),fmt,s,tag); f_out=fopen(out_filename,"wt"); free(out_filename); printf("Constructing percentiles for forecasts - regime %d\n",s); @@ -469,7 +471,7 @@ int main(int nargs, char **args) else { fmt="forecasts_percentiles_%s.prn"; - sprintf(out_filename=(char*)malloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); + sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); f_out=fopen(out_filename,"wt"); free(out_filename); printf("Constructing percentiles for forecasts - %d draws of shocks/regimes per posterior value\n",draws); 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 79ef5bd4a..70241c927 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -6,6 +6,8 @@ #include "dw_histogram.h" #include "dw_error.h" +#include "modify_for_mex.h" + static void Resize(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); static void AddObservationVariable(PRECISION x, int *h, PRECISION *min, PRECISION *max, int intervals); static void AddObservationFixed(PRECISION x, int *low, int *h, int *high, PRECISION min, PRECISION max, int intervals); @@ -37,23 +39,23 @@ TMatrixHistogram *CreateMatrixHistogram(int rows, int cols, int intervals, int t int i, j; TMatrixHistogram *h; - if (!(h=(TMatrixHistogram *)malloc(sizeof(TMatrixHistogram)))) dw_Error(MEM_ERR); + if (!(h=(TMatrixHistogram *)swzMalloc(sizeof(TMatrixHistogram)))) dw_Error(MEM_ERR); - if (!(h->freq=(int***)malloc(rows*sizeof(int**)))) dw_Error(MEM_ERR); + if (!(h->freq=(int***)swzMalloc(rows*sizeof(int**)))) dw_Error(MEM_ERR); for (i=rows-1; i >= 0; i--) { - if (!(h->freq[i]=(int**)malloc(cols*sizeof(int*)))) dw_Error(MEM_ERR); + if (!(h->freq[i]=(int**)swzMalloc(cols*sizeof(int*)))) dw_Error(MEM_ERR); for (j=cols-1; j >= 0; j--) - if (!(h->freq[i][j]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->freq[i][j]=(int*)swzMalloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); } - if (!(h->low=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + if (!(h->low=(int**)swzMalloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); for (i=rows-1; i >= 0; i--) - if (!(h->low[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->low[i]=(int*)swzMalloc(cols*sizeof(int)))) dw_Error(MEM_ERR); - if (!(h->high=(int**)malloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); + if (!(h->high=(int**)swzMalloc(rows*sizeof(int*)))) dw_Error(MEM_ERR); for (i=rows-1; i >= 0; i--) - if (!(h->high[i]=(int*)malloc(cols*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->high[i]=(int*)swzMalloc(cols*sizeof(int)))) dw_Error(MEM_ERR); h->Min=CreateMatrix(rows,cols); h->Max=CreateMatrix(rows,cols); @@ -193,15 +195,15 @@ TVectorHistogram *CreateVectorHistogram(int dim, int intervals, int type) int i; TVectorHistogram *h; - if (!(h=(TVectorHistogram *)malloc(sizeof(TVectorHistogram)))) + if (!(h=(TVectorHistogram *)swzMalloc(sizeof(TVectorHistogram)))) dw_Error(MEM_ERR); - if (!(h->freq=(int**)malloc(dim*sizeof(int*)))) dw_Error(MEM_ERR); + if (!(h->freq=(int**)swzMalloc(dim*sizeof(int*)))) dw_Error(MEM_ERR); for (i=dim-1; i >= 0; i--) - if (!(h->freq[i]=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->freq[i]=(int*)swzMalloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); - if (!(h->low=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); - if (!(h->high=(int*)malloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->low=(int*)swzMalloc(dim*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->high=(int*)swzMalloc(dim*sizeof(int)))) dw_Error(MEM_ERR); h->Min=CreateVector(dim); h->Max=CreateVector(dim); @@ -333,9 +335,9 @@ TScalarHistogram *CreateScalarHistogram(int intervals, int type) { TScalarHistogram *h; - if (!(h=(TScalarHistogram *)malloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR); + if (!(h=(TScalarHistogram *)swzMalloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR); - if (!(h->freq=(int*)malloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); + if (!(h->freq=(int*)swzMalloc(intervals*sizeof(int)))) dw_Error(MEM_ERR); h->intervals=intervals; h->sample_size=0; From 820aca66839d755eefa2a5eb107095051576a6ff Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 18:01:47 +0200 Subject: [PATCH 06/12] SWZ: replace free with swzFree --- matlab/swz/c-code/sbvar/var/forecast.c | 18 +++++------ .../DWCcode/histogram/dw_histogram.c | 30 +++++++++---------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index d98baebd0..b1b9f6ddd 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -119,7 +119,7 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste ERROR_EXIT: FreeMatrixHistogram(histogram); FreeMatrix(forecast); - free(S); + swzFree(S); FreeVector(initial); FreeVector(prob); FreeVector(init_prob); @@ -229,7 +229,7 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, ERROR_EXIT: FreeMatrixHistogram(histogram); FreeMatrix(forecast); - free(S); + swzFree(S); FreeVector(initial); FreeVector(prob); FreeVector(init_prob); @@ -362,9 +362,9 @@ int main(int nargs, char **args) Read_VAR_Parameters((FILE*)NULL,parm,head,model); p=(T_VAR_Parameters*)(model->theta); - free(spec); - free(head); - free(parm); + swzFree(spec); + swzFree(head); + swzFree(parm); //============================= Compute forecasts ============================= @@ -374,7 +374,7 @@ int main(int nargs, char **args) /* fmt="forecasts_mean_%s.prn"; */ /* sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); */ /* f_out=fopen(out_filename,"wt"); */ - /* free(out_filename); */ + /* 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++) */ @@ -452,7 +452,7 @@ int main(int nargs, char **args) 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"); - free(out_filename); + 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); @@ -463,7 +463,7 @@ int main(int nargs, char **args) 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"); - free(out_filename); + 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); @@ -473,7 +473,7 @@ int main(int nargs, char **args) fmt="forecasts_percentiles_%s.prn"; sprintf(out_filename=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); f_out=fopen(out_filename,"wt"); - free(out_filename); + 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); 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 70241c927..627a18b21 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c @@ -81,17 +81,17 @@ void FreeMatrixHistogram(TMatrixHistogram *h) int i, j; for (i=h->rows-1; i >= 0; i--) { - for (j=h->cols-1; j >= 0; j--) free(h->freq[i][j]); - free(h->freq[i]); + for (j=h->cols-1; j >= 0; j--) swzFree(h->freq[i][j]); + swzFree(h->freq[i]); } - free(h->freq); - for (i=h->rows-1; i >= 0; i--) free(h->low[i]); - free(h->low); - for (i=h->rows-1; i >= 0; i--) free(h->high[i]); - free(h->high); + swzFree(h->freq); + for (i=h->rows-1; i >= 0; i--) swzFree(h->low[i]); + swzFree(h->low); + for (i=h->rows-1; i >= 0; i--) swzFree(h->high[i]); + swzFree(h->high); FreeMatrix(h->Min); FreeMatrix(h->Max); - free(h); + swzFree(h); } void AddMatrixObservation(TMatrix X, TMatrixHistogram *h) @@ -226,13 +226,13 @@ void SetMaxMinVectorHistogram(TVector Min, TVector Max, TVectorHistogram *h) void FreeVectorHistogram(TVectorHistogram *h) { int i; - for (i=h->dim-1; i >= 0; i--) free(h->freq[i]); - free(h->freq); - free(h->low); - free(h->high); + for (i=h->dim-1; i >= 0; i--) swzFree(h->freq[i]); + swzFree(h->freq); + swzFree(h->low); + swzFree(h->high); FreeVector(h->Min); FreeVector(h->Max); - free(h); + swzFree(h); } void AddVectorObservation(TVector x, TVectorHistogram *h) @@ -355,8 +355,8 @@ void SetMaxMinScalarHistogram(PRECISION Min, PRECISION Max, TScalarHistogram *h) void FreeScalarHistogram(TScalarHistogram *h) { - free(h->freq); - free(h); + swzFree(h->freq); + swzFree(h); } void AddScalarObservation(PRECISION x, TScalarHistogram *h) From 94217de476df2ab76674cebf6d4aae39a9160dfe Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 18:03:57 +0200 Subject: [PATCH 07/12] SWZ: replace exit with swzExit --- matlab/swz/c-code/sbvar/var/forecast.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index b1b9f6ddd..5dfdd4000 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -345,7 +345,7 @@ int main(int nargs, char **args) " -fh : parameter header (Posterior mode: )\n" " -horizon : horizon for the forecast (12)\n" ); - exit(1); + swzExit(1); } if (!parm) @@ -392,7 +392,7 @@ int main(int nargs, char **args) if (!(posterior_file=fopen(post,"rt"))) { printf("Unable to open draws file: %s\n",post); - exit(0); + swzExit(0); } // Get thinning factor from command line @@ -436,13 +436,13 @@ int main(int nargs, char **args) { FreeVector(percentiles); printf("forecasting command line: Error parsing percentiles\n"); - exit(0); + swzExit(0); } } else { printf("forecasting command line(): Error parsing percentiles\n"); - exit(0); + swzExit(0); } if (dw_FindArgument_String(nargs,args,"regimes") != -1) From 8cf788af44f5128e0b58b6d5ef75a558b76dcbd3 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Aug 2010 18:08:36 +0200 Subject: [PATCH 08/12] SWZ: replace fprintf(stderr, with swz_fprintf_err( --- matlab/swz/c-code/sbvar/var/forecast.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index 5dfdd4000..0d7d90e3d 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -337,8 +337,8 @@ int main(int nargs, char **args) if (!spec) { - fprintf(stderr,"No specification filename given\n"); - fprintf(stderr,"Command line syntax:\n" + 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" From 8b782e7b005e349ff7ef62570c2da17668171ae2 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 30 Aug 2010 18:20:33 +0200 Subject: [PATCH 09/12] SWZ: replace matrix.h with swzmatrix.h --- matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h index 64e24a2e3..619b5265a 100644 --- a/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h +++ b/matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h @@ -1,7 +1,7 @@ #ifndef __HISTOGRAMS__ #define __HISTOGRAMS__ -#include "matrix.h" +#include "swzmatrix.h" #define HISTOGRAM_FIXED 1 #define HISTOGRAM_VARIABLE 2 From 60a3c2cad478296d024dcc4e07f3528c06acb137 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 30 Aug 2010 18:56:26 +0200 Subject: [PATCH 10/12] SWZ: replace tabs with spaces --- matlab/swz/c-code/sbvar/switching/switchio.c | 64 ++--- matlab/swz/c-code/sbvar/var/VARbase.c | 108 ++++----- matlab/swz/c-code/sbvar/var/VARio.c | 32 +-- matlab/swz/c-code/sbvar/var/forecast.c | 220 +++++++++--------- .../DWCcode/histogram/dw_histogram.c | 114 ++++----- 5 files changed, 269 insertions(+), 269 deletions(-) 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; + } } } From 137c4cf6b112ee08af102b9a4c55159f18d33752 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 30 Aug 2010 18:57:47 +0200 Subject: [PATCH 11/12] SWZ: make comments conform to ansi C standard --- matlab/swz/c-code/sbvar/switching/switchio.c | 10 +-- matlab/swz/c-code/sbvar/var/VARbase.c | 8 +-- matlab/swz/c-code/sbvar/var/VARbase.h | 4 +- matlab/swz/c-code/sbvar/var/VARio.c | 12 ++-- matlab/swz/c-code/sbvar/var/forecast.c | 74 ++++++++++---------- 5 files changed, 54 insertions(+), 54 deletions(-) diff --git a/matlab/swz/c-code/sbvar/switching/switchio.c b/matlab/swz/c-code/sbvar/switching/switchio.c index 7f246f45d..fe1efd679 100644 --- a/matlab/swz/c-code/sbvar/switching/switchio.c +++ b/matlab/swz/c-code/sbvar/switching/switchio.c @@ -733,7 +733,7 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) } else { - // Read transition matrix +/* // Read transition matrix ansi-c*/ Q=CreateMatrix(sv->nbasestates,sv->nbasestates); for (j=0; j < ColM(Q); j++) for (i=0; i < RowM(Q); i++) @@ -743,7 +743,7 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) return 0; } - // Scale the columns of Q - loose requirement on sumation to one +/* // Scale the columns of Q - loose requirement on sumation to one ansi-c*/ for (j=sv->nbasestates-1; j >= 0; j--) { for (sum=0.0, i=sv->nbasestates-1; i >= 0; i--) @@ -765,13 +765,13 @@ int ReadBaseTransitionMatricesFlat_SV(FILE *f_in, TMarkovStateVariable *sv) ElementM(Q,i,j)*=sum; } - // Convert base transition matrix to full transition matrix. +/* // Convert base transition matrix to full transition matrix. ansi-c*/ ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded); - // Free Q +/* // Free Q ansi-c*/ FreeMatrix(Q); - // Update +/* // Update ansi-c*/ if (!Update_B_from_Q_SV(sv)) { dw_UserError("Transition matrices do not satisfy restrictions"); diff --git a/matlab/swz/c-code/sbvar/var/VARbase.c b/matlab/swz/c-code/sbvar/var/VARbase.c index 186ef20c3..aa8414542 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.c +++ b/matlab/swz/c-code/sbvar/var/VARbase.c @@ -2418,17 +2418,17 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s TVector x, y; int i, t; - // allocate forecast if necessary +/* // allocate forecast if necessary ansi-c*/ if (!forecast && !(forecast=CreateMatrix(horizon,p->nvars))) return (TMatrix)NULL; - // allocate memory +/* // allocate memory ansi-c*/ y=CreateVector(p->nvars); x=EquateVector((TVector)NULL,initial); A0=MakeA0_All((TMatrix*)NULL,p); Aplus=MakeAplus_All((TMatrix*)NULL,p); - // forecast +/* // forecast ansi-c*/ for (t=0; t < horizon; t++) { ProductVM(y,x,Aplus[S[t]]); @@ -2445,7 +2445,7 @@ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *s memcpy(pElementV(x),pElementV(y),p->nvars*sizeof(PRECISION)); } - // free memory +/* // free memory ansi-c*/ dw_FreeArray(Aplus); dw_FreeArray(A0); FreeVector(x); diff --git a/matlab/swz/c-code/sbvar/var/VARbase.h b/matlab/swz/c-code/sbvar/var/VARbase.h index e221daa1c..a4f2d28c2 100644 --- a/matlab/swz/c-code/sbvar/var/VARbase.h +++ b/matlab/swz/c-code/sbvar/var/VARbase.h @@ -201,8 +201,8 @@ void Draw_lambda(TStateModel *model); /* Forecasts */ TMatrix forecast_base(TMatrix forecast, int horizon, TVector initial, TVector *shocks, int *S, TStateModel *model); -//TVector* mean_conditional_forecast(TVector *F, PRECISION ***y, int h, int t0, TStateModel *model); -//TVector* mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model); +/* //TVector* mean_conditional_forecast(TVector *F, PRECISION ***y, int h, int t0, TStateModel *model); ansi-c*/ +/* //TVector* mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model); ansi-c*/ /* Utilities */ void ComputeDotProducts_All(T_VAR_Parameters *p); diff --git a/matlab/swz/c-code/sbvar/var/VARio.c b/matlab/swz/c-code/sbvar/var/VARio.c index f60c370c3..9fb0528ed 100644 --- a/matlab/swz/c-code/sbvar/var/VARio.c +++ b/matlab/swz/c-code/sbvar/var/VARio.c @@ -570,12 +570,12 @@ int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model) int i, j, s, rtrn=0; T_VAR_Parameters *p=(T_VAR_Parameters*)(model->theta); - // Allocate memory +/* // Allocate memory ansi-c*/ A0=dw_CreateArray_matrix(p->nstates); Aplus=dw_CreateArray_matrix(p->nstates); Zeta=dw_CreateArray_vector(p->nstates); - // Read File +/* // Read File ansi-c*/ for (s=0; s < p->nstates; s++) { A0[s]=CreateMatrix(p->nvars,p->nvars); @@ -599,7 +599,7 @@ int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model) goto ERROR; } - // Set A0, Aplus, and Zeta +/* // Set A0, Aplus, and Zeta ansi-c*/ for (j=0; j < p->nvars; j++) for (s=0; s < p->nstates; s++) { @@ -612,18 +612,18 @@ int Read_VAR_ParametersFlat(FILE *f_in, TStateModel *model) p->Zeta[j][p->var_states[j][s]]=ElementV(Zeta[s],j); } - // Update b0, bplus, lambda, psi +/* // Update b0, bplus, lambda, psi ansi-c*/ Update_b0_bplus_from_A0_Aplus(p); if ((p->Specification & SPEC_SIMS_ZHA) == SPEC_SIMS_ZHA) Update_lambda_psi_from_bplus(p); - // Flags and notification that the VAR parameters have changed +/* // Flags and notification that the VAR parameters have changed ansi-c*/ p->valid_parameters=1; ThetaChanged(model); rtrn=1; ERROR: - // Free memory +/* // Free memory ansi-c*/ dw_FreeArray(A0); dw_FreeArray(Aplus); dw_FreeArray(Zeta); diff --git a/matlab/swz/c-code/sbvar/var/forecast.c b/matlab/swz/c-code/sbvar/var/forecast.c index 5f519e75b..1493a053d 100644 --- a/matlab/swz/c-code/sbvar/var/forecast.c +++ b/matlab/swz/c-code/sbvar/var/forecast.c @@ -37,14 +37,14 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste TMatrix forecast; TMatrixHistogram *histogram; - // quick check of passed parameters +/* // quick check of passed parameters ansi-c*/ if (!f_out || !percentiles || (draws <= 0) || (T < 0) || (h < 0) || !model) return 0; p=(T_VAR_Parameters*)(model->theta); if (T > p->nobs) return 0; - // allocate memory +/* // allocate memory ansi-c*/ S=(int*)swzMalloc(h*sizeof(int)); forecast=CreateMatrix(h,p->nvars); histogram=CreateMatrixHistogram(h,p->nvars,100,HISTOGRAM_VARIABLE); @@ -54,13 +54,13 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste init_prob=CreateVector(p->nstates); prob=CreateVector(p->nstates); - // Initial value +/* // Initial value ansi-c*/ EquateVector(initial,p->X[T]); i=0; while (!done) { - // Read parameters and push them into model +/* // Read parameters and push them into model ansi-c*/ if (!posterior_file) done=1; else @@ -78,30 +78,30 @@ int forecast_percentile(FILE *f_out, TVector percentiles, int draws, FILE *poste if (done != 2) { - // Get filtered probability at time T +/* // Get filtered probability at time T ansi-c*/ 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 +/* // Draw time T regime ansi-c*/ m=DrawDiscrete(init_prob); - // Draw regimes from time T+1 through T+h inclusive +/* // Draw regimes from time T+1 through T+h inclusive ansi-c*/ 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 ansi-c*/ + for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); /* InitializeVector(shocks[i],0.0); ansi-c*/ - // Compute forecast +/* // Compute forecast ansi-c*/ if (!forecast_base(forecast,h,initial,shocks,S,model)) goto ERROR_EXIT; - // Accumulate impulse response +/* // Accumulate impulse response ansi-c*/ AddMatrixObservation(forecast,histogram); } } @@ -160,14 +160,14 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, TMatrix forecast; TMatrixHistogram *histogram; - // quick check of passed parameters +/* // quick check of passed parameters ansi-c*/ if (!f_out || !percentiles || (draws <= 0) || (T < 0) || (h < 0) || !model) return 0; p=(T_VAR_Parameters*)(model->theta); if (T > p->nobs) return 0; - // allocate memory +/* // allocate memory ansi-c*/ S=(int*)swzMalloc(h*sizeof(int)); for (i=0; i < h; i++) S[i]=s; forecast=CreateMatrix(h,p->nvars); @@ -178,13 +178,13 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, init_prob=CreateVector(p->nstates); prob=CreateVector(p->nstates); - // Initial value +/* // Initial value ansi-c*/ EquateVector(initial,p->X[T]); i=0; while (!done) { - // Read parameters and push them into model +/* // Read parameters and push them into model ansi-c*/ if (!posterior_file) done=1; else @@ -204,14 +204,14 @@ int forecast_percentile_regime(FILE *f_out, TVector percentiles, int draws, { for (k=draws; k > 0; k--) { - // Draw shocks - for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); // InitializeVector(shocks[i],0.0); +/* // Draw shocks ansi-c*/ + for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); /* InitializeVector(shocks[i],0.0); ansi-c*/ - // Compute forecast +/* // Compute forecast ansi-c*/ if (!forecast_base(forecast,h,initial,shocks,S,model)) goto ERROR_EXIT; - // Accumulate impulse response +/* // Accumulate impulse response ansi-c*/ AddMatrixObservation(forecast,histogram); } } @@ -306,33 +306,33 @@ int main(int nargs, char **args) int s, horizon, thin, draws, i, j, n; FILE *f_out, *posterior_file; - // specification filename +/* // specification filename ansi-c*/ if (buffer=dw_ParseString_String(nargs,args,"fs",(char*)NULL)) strcpy(spec=(char*)swzMalloc(strlen(buffer)+1),buffer); - // parameter filename +/* // parameter filename ansi-c*/ if (buffer=dw_ParseString_String(nargs,args,"fp",(char*)NULL)) strcpy(parm=(char*)swzMalloc(strlen(buffer)+1),buffer); - // header +/* // header ansi-c*/ if (buffer=dw_ParseString_String(nargs,args,"ph",(char*)NULL)) strcpy(head=(char*)swzMalloc(strlen(buffer)+1),buffer); - // file tag +/* // file tag ansi-c*/ if (tag=dw_ParseString_String(nargs,args,"ft",(char*)NULL)) { fmt="est_final_%s.dat"; - // specification filename +/* // specification filename ansi-c*/ if (!spec) sprintf(spec=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); - // parameter filename +/* // parameter filename ansi-c*/ if (!parm) sprintf(parm=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); } - // horizon +/* // horizon ansi-c*/ horizon=dw_ParseInteger_String(nargs,args,"horizon",12); if (!spec) @@ -366,9 +366,9 @@ int main(int nargs, char **args) swzFree(head); swzFree(parm); - //============================= Compute forecasts ============================= +/* //============================= Compute forecasts ============================= ansi-c*/ - // Mean forecast +/* // Mean forecast ansi-c*/ /* if (dw_FindArgument_String(nargs,args,"mean") != -1) */ /* { */ /* fmt="forecasts_mean_%s.prn"; */ @@ -383,10 +383,10 @@ int main(int nargs, char **args) /* return; */ /* } */ - // Parameter uncertainty +/* // Parameter uncertainty ansi-c*/ if (dw_FindArgument_String(nargs,args,"parameter_uncertainity") != -1) { - // Open posterior draws file +/* // Open posterior draws file ansi-c*/ fmt="draws_%s.dat"; sprintf(post=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag); if (!(posterior_file=fopen(post,"rt"))) @@ -395,25 +395,25 @@ int main(int nargs, char **args) swzExit(0); } - // Get thinning factor from command line +/* // Get thinning factor from command line ansi-c*/ thin=dw_ParseInteger_String(nargs,args,"thin",1); - // Get shocks_per_parameter from command line +/* // Get shocks_per_parameter from command line ansi-c*/ draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",1); } else { - // Using posterior estimate +/* // Using posterior estimate ansi-c*/ posterior_file=(FILE*)NULL; - // thinning factor not used +/* // thinning factor not used ansi-c*/ thin=1; - // Get shocks_per_parameter from command line +/* // Get shocks_per_parameter from command line ansi-c*/ draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",10000); } - // Setup percentiles +/* // Setup percentiles ansi-c*/ if ((i=dw_FindArgument_String(nargs,args,"percentiles")) == -1) if (dw_FindArgument_String(nargs,args,"error_bands") == -1) { @@ -481,7 +481,7 @@ int main(int nargs, char **args) if (posterior_file) fclose(posterior_file); FreeVector(percentiles); - //============================================================================= +/* //============================================================================= ansi-c*/ return 0; } From 18c947917617923b5b5150927c86512555d29f6b Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 30 Aug 2010 18:47:32 +0200 Subject: [PATCH 12/12] SWZ: include in build system --- matlab/swz/c-code/Makefile | 13 ++++++++----- mex/build/swz.am | 10 ++++++---- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/matlab/swz/c-code/Makefile b/matlab/swz/c-code/Makefile index e17d93406..c99c196e8 100644 --- a/matlab/swz/c-code/Makefile +++ b/matlab/swz/c-code/Makefile @@ -191,6 +191,7 @@ MATRIX_DIR = $(WORK_DIR)/utilities/DWCcode/matrix ERROR_DIR = $(WORK_DIR)/utilities/DWCcode/error ARRAY_DIR = $(WORK_DIR)/utilities/DWCcode/arrays ASCII_DIR = $(WORK_DIR)/utilities/DWCcode/ascii +HIST_DIR = $(WORK_DIR)/utilities/DWCcode/histogram STAT_DIR = $(WORK_DIR)/utilities/DWCcode/stat SPHERICAL_DIR = $(WORK_DIR)/utilities/DWCcode/spherical SORT_DIR = $(WORK_DIR)/utilities/DWCcode/sort @@ -199,9 +200,9 @@ VAR_DIR = $(WORK_DIR)/sbvar/var ################################################################################# # DW FILES -INCLUDE_DIR := -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(STAT_DIR) -I$(SPHERICAL_DIR) -I$(SORT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR) $(INCLUDE_DIR) -VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(STAT_DIR) $(SPHERICAL_DIR) $(SORT_DIR) $(SWITCH_DIR) $(VAR_DIR) -OBJS := $(OBJS) bmatrix.o swzmatrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_matrix_sort.o dw_ascii.o dw_parse_cmd.o +INCLUDE_DIR := -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(HIST_DIR) -I$(STAT_DIR) -I$(SPHERICAL_DIR) -I$(SORT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR) $(INCLUDE_DIR) +VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(HIST_DIR) $(STAT_DIR) $(SPHERICAL_DIR) $(SORT_DIR) $(SWITCH_DIR) $(VAR_DIR) +OBJS := $(OBJS) bmatrix.o swzmatrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_matrix_sort.o dw_ascii.o dw_parse_cmd.o dw_histogram.o # MEX INCLUDE_DIR := -I$(WORK_DIR)/mex $(INCLUDE_DIR) @@ -221,13 +222,13 @@ OBJS_INIT := $(OBJS) create_init_file.o switch.o switchio.o VARbase.o VARio.o VA OBJS_MHM_1 := $(OBJS) mhm_VAR_main_1.o mhm_VAR.o VARbase.o VARio.o command_line_VAR.o switch.o switchio.o OBJS_MHM_2 := $(OBJS) mhm_VAR_main_2.o spherical.o VARbase.o VARio.o switch.o switchio.o mhm_VAR.o OBJS_PROBA := $(OBJS) probabilities.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o - +OBJS_FORECAST := $(OBJS) forecast.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o ################################################################################# # OUTPUT -all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities +all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities $(OUT_DIR)/sbvar_forecast $(OUT_DIR)/sbvar_draws: $(OBJS_DRAWS) @@ -248,6 +249,8 @@ $(OUT_DIR)/sbvar_mhm_2: $(OBJS_MHM_2) $(OUT_DIR)/sbvar_probabilities: $(OBJS_PROBA) $(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_probabilities +$(OUT_DIR)/sbvar_forecast: $(OBJS_FORECAST) + $(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_forecast %.o : %.c diff --git a/mex/build/swz.am b/mex/build/swz.am index 80925a984..83c157698 100644 --- a/mex/build/swz.am +++ b/mex/build/swz.am @@ -1,18 +1,18 @@ SWZ_SRC_BASEDIR = ../../../../matlab/swz/c-code -SWZ_SRC_DIRS = $(SWZ_SRC_BASEDIR)/utilities/TZCcode $(SWZ_SRC_BASEDIR)/utilities/DWCcode/matrix $(SWZ_SRC_BASEDIR)/utilities/DWCcode/error $(SWZ_SRC_BASEDIR)/utilities/DWCcode/arrays $(SWZ_SRC_BASEDIR)/utilities/DWCcode/ascii $(SWZ_SRC_BASEDIR)/utilities/DWCcode/stat $(SWZ_SRC_BASEDIR)/utilities/DWCcode/spherical $(SWZ_SRC_BASEDIR)/utilities/DWCcode/sort $(SWZ_SRC_BASEDIR)/sbvar/switching $(SWZ_SRC_BASEDIR)/sbvar/var $(SWZ_SRC_BASEDIR)/mex +SWZ_SRC_DIRS = $(SWZ_SRC_BASEDIR)/utilities/TZCcode $(SWZ_SRC_BASEDIR)/utilities/DWCcode/matrix $(SWZ_SRC_BASEDIR)/utilities/DWCcode/error $(SWZ_SRC_BASEDIR)/utilities/DWCcode/arrays $(SWZ_SRC_BASEDIR)/utilities/DWCcode/ascii $(SWZ_SRC_BASEDIR)/utilities/DWCcode/histogram $(SWZ_SRC_BASEDIR)/utilities/DWCcode/stat $(SWZ_SRC_BASEDIR)/utilities/DWCcode/spherical $(SWZ_SRC_BASEDIR)/utilities/DWCcode/sort $(SWZ_SRC_BASEDIR)/sbvar/switching $(SWZ_SRC_BASEDIR)/sbvar/var $(SWZ_SRC_BASEDIR)/mex vpath %.c $(SWZ_SRC_DIRS) -CPPFLAGS += -DINTELCMATHLIBRARY $(GSL_CPPFLAGS) -I$(SWZ_SRC_BASEDIR)/utilities/TZCcode -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/matrix -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/error -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/arrays -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/ascii -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/stat -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/spherical -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/sort -I$(SWZ_SRC_BASEDIR)/sbvar/switching -I$(SWZ_SRC_BASEDIR)/sbvar/var -I$(SWZ_SRC_BASEDIR)/mex +CPPFLAGS += -DINTELCMATHLIBRARY $(GSL_CPPFLAGS) -I$(SWZ_SRC_BASEDIR)/utilities/TZCcode -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/matrix -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/error -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/arrays -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/ascii -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/histogram -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/stat -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/spherical -I$(SWZ_SRC_BASEDIR)/utilities/DWCcode/sort -I$(SWZ_SRC_BASEDIR)/sbvar/switching -I$(SWZ_SRC_BASEDIR)/sbvar/var -I$(SWZ_SRC_BASEDIR)/mex LIBS += $(GSL_LIBS) LDFLAGS += $(GSL_LDFLAGS) -noinst_PROGRAMS = mex_sbvar_init_file mex_sbvar_estimation mex_sbvar_mhm_1 mex_sbvar_mhm_2 mex_sbvar_probabilities mex_sbvar_draws +noinst_PROGRAMS = mex_sbvar_init_file mex_sbvar_estimation mex_sbvar_mhm_1 mex_sbvar_mhm_2 mex_sbvar_probabilities mex_sbvar_draws mex_sbvar_forecast common_mex = mex_top_level.c modify_for_mex.c -swz_common = bmatrix.c swzmatrix.c dw_error.c dw_rand.c dw_matrix_rand.c dw_array.c dw_matrix_array.c dw_matrix_sort.c dw_ascii.c dw_parse_cmd.c +swz_common = bmatrix.c swzmatrix.c dw_error.c dw_rand.c dw_matrix_rand.c dw_array.c dw_matrix_array.c dw_matrix_sort.c dw_ascii.c dw_parse_cmd.c dw_histogram.c swz_tao = tzmatlab.c mathlib.c cstz_dw.c nodist_mex_sbvar_init_file_SOURCES = $(common_mex) $(swz_common) create_init_file.c switch.c switchio.c VARbase.c VARio.c VARio_matlab.c @@ -21,3 +21,5 @@ nodist_mex_sbvar_mhm_1_SOURCES = $(common_mex) $(swz_common) mhm_VAR_mai nodist_mex_sbvar_mhm_2_SOURCES = $(common_mex) $(swz_common) mhm_VAR_main_2.c spherical.c VARbase.c VARio.c switch.c switchio.c mhm_VAR.c nodist_mex_sbvar_probabilities_SOURCES = $(common_mex) $(swz_common) probabilities.c switch.c switchio.c VARbase.c VARio.c command_line_VAR.c nodist_mex_sbvar_draws_SOURCES = $(common_mex) $(swz_common) PrintDraws.c switch.c switchio.c VARbase.c VARio.c command_line_VAR.c +nodist_mex_sbvar_forecast_SOURCES = $(common_mex) $(swz_common) forecast.c switch.c switchio.c VARbase.c VARio.c command_line_VAR.c +