stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of kirikou.dynare.org:/srv/d_kirikou/git/dynare

time-shift
Stéphane Adjemian (Scylla) 2010-09-01 22:14:32 +02:00
parent 6870010e9a 18c9479176
commit ded34f990b
11 changed files with 1773 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
matlab/swz/c-code/Makefile
View File

@ -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

77
matlab/swz/c-code/sbvar/switching/switchio.c
View File

@ -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 ansi-c*/
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 ansi-c*/
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. ansi-c*/
ConvertBaseTransitionMatrix(sv->Q,Q,sv->nlags_encoded);
/* // Free Q ansi-c*/
FreeMatrix(Q);
/* // Update ansi-c*/
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;

2
matlab/swz/c-code/sbvar/switching/switchio.h
View File

@ -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);
/*

286
matlab/swz/c-code/sbvar/var/VARbase.c
View File

@ -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 ansi-c*/
if (!forecast && !(forecast=CreateMatrix(horizon,p->nvars)))
return (TMatrix)NULL;
/* // 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 ansi-c*/
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 ansi-c*/
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

9
matlab/swz/c-code/sbvar/var/VARbase.h
View File

@ -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); ansi-c*/
/* //TVector* mean_unconditional_forecast(TVector *F, int h, int t0, TStateModel *model); ansi-c*/
/* 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);

68
matlab/swz/c-code/sbvar/var/VARio.c
View File

@ -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 ansi-c*/
A0=dw_CreateArray_matrix(p->nstates);
Aplus=dw_CreateArray_matrix(p->nstates);
Zeta=dw_CreateArray_vector(p->nstates);
/* // Read File ansi-c*/
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 ansi-c*/
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 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 ansi-c*/
p->valid_parameters=1;
ThetaChanged(model);
rtrn=1;
ERROR:
/* // Free memory ansi-c*/
dw_FreeArray(A0);
dw_FreeArray(Aplus);
dw_FreeArray(Zeta);
return rtrn;
}
/*
For each state the VAR parameters are printed as follows
A0 (by columns)

1
matlab/swz/c-code/sbvar/var/VARio.h
View File

@ -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);

488
matlab/swz/c-code/sbvar/var/forecast.c Normal file
View File

@ -0,0 +1,488 @@
#include "switch.h"
#include "switchio.h"
#include "VARio.h"
#include "dw_parse_cmd.h"
#include "dw_ascii.h"
#include "dw_histogram.h"
#include <stdlib.h>
#include <string.h>
#include "modify_for_mex.h"
/*
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 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 ansi-c*/
S=(int*)swzMalloc(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 ansi-c*/
EquateVector(initial,p->X[T]);
i=0;
while (!done)
{
/* // Read parameters and push them into model ansi-c*/
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 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 ansi-c*/
m=DrawDiscrete(init_prob);
/* // 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 ansi-c*/
for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); /* InitializeVector(shocks[i],0.0); ansi-c*/
/* // Compute forecast ansi-c*/
if (!forecast_base(forecast,h,initial,shocks,S,model))
goto ERROR_EXIT;
/* // Accumulate impulse response ansi-c*/
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);
swzFree(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 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 ansi-c*/
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);
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 ansi-c*/
EquateVector(initial,p->X[T]);
i=0;
while (!done)
{
/* // Read parameters and push them into model ansi-c*/
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 ansi-c*/
for (j=h-1; j >= 0; j--) dw_NormalVector(shocks[j]); /* InitializeVector(shocks[i],0.0); ansi-c*/
/* // Compute forecast ansi-c*/
if (!forecast_base(forecast,h,initial,shocks,S,model))
goto ERROR_EXIT;
/* // Accumulate impulse response ansi-c*/
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);
swzFree(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 <filename tag>
If this argument exists, then the following is attempted:
specification file name = est_final_<tag>.dat
output file name = ir_<tag>_regime_<k>.dat
parameters file name = est_final_<tag>.dat
header = "Posterior mode: "
-fs <filename>
If this argument exists, then the specification file name is <filename>.
The argument -fs takes precedence over -ft.
-fp <filename>
If this argument exists, then the parameters file name is <filename>. The
argument -fp takes precedence over -ft. The default value is the filename
associated with the argument -fs.
-ph <header>
If this argument exists, then the header for the parameters file is
<header>. The default value is "Posterior mode: ".
-horizon <integer>
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 <integer>
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 <integer>
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 ansi-c*/
if (buffer=dw_ParseString_String(nargs,args,"fs",(char*)NULL))
strcpy(spec=(char*)swzMalloc(strlen(buffer)+1),buffer);
/* // parameter filename ansi-c*/
if (buffer=dw_ParseString_String(nargs,args,"fp",(char*)NULL))
strcpy(parm=(char*)swzMalloc(strlen(buffer)+1),buffer);
/* // header ansi-c*/
if (buffer=dw_ParseString_String(nargs,args,"ph",(char*)NULL))
strcpy(head=(char*)swzMalloc(strlen(buffer)+1),buffer);
/* // file tag ansi-c*/
if (tag=dw_ParseString_String(nargs,args,"ft",(char*)NULL))
{
fmt="est_final_%s.dat";
/* // specification filename ansi-c*/
if (!spec)
sprintf(spec=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag);
/* // parameter filename ansi-c*/
if (!parm)
sprintf(parm=(char*)swzMalloc(strlen(fmt) + strlen(tag) - 1),fmt,tag);
}
/* // horizon ansi-c*/
horizon=dw_ParseInteger_String(nargs,args,"horizon",12);
if (!spec)
{
swz_fprintf_err("No specification filename given\n");
swz_fprintf_err("Command line syntax:\n"
" -ft : file tag\n"
" -fs : specification filename (est_final_<tag>.dat)\n"
" -fp : parameters filename (specification filename)\n"
" -fh : parameter header (Posterior mode: )\n"
" -horizon : horizon for the forecast (12)\n"
);
swzExit(1);
}
if (!parm)
strcpy(parm=(char*)swzMalloc(strlen(spec)+1),spec);
if (!head)
{
buffer="Posterior mode: ";
strcpy(head=(char*)swzMalloc(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);
swzFree(spec);
swzFree(head);
swzFree(parm);
/* //============================= Compute forecasts ============================= ansi-c*/
/* // Mean forecast ansi-c*/
/* if (dw_FindArgument_String(nargs,args,"mean") != -1) */
/* { */
/* fmt="forecasts_mean_%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 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 ansi-c*/
if (dw_FindArgument_String(nargs,args,"parameter_uncertainity") != -1)
{
/* // 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")))
{
printf("Unable to open draws file: %s\n",post);
swzExit(0);
}
/* // Get thinning factor from command line ansi-c*/
thin=dw_ParseInteger_String(nargs,args,"thin",1);
/* // Get shocks_per_parameter from command line ansi-c*/
draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",1);
}
else
{
/* // Using posterior estimate ansi-c*/
posterior_file=(FILE*)NULL;
/* // thinning factor not used ansi-c*/
thin=1;
/* // Get shocks_per_parameter from command line ansi-c*/
draws=dw_ParseInteger_String(nargs,args,"shocks_per_parameter",10000);
}
/* // Setup percentiles ansi-c*/
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");
swzExit(0);
}
}
else
{
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);
}
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);
}
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);
}
if (posterior_file) fclose(posterior_file);
FreeVector(percentiles);
/* //============================================================================= ansi-c*/
return 0;
}

751
matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.c Normal file
View File

@ -0,0 +1,751 @@
#include <math.h>
#include <string.h>
#include <stdlib.h>
#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);
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 *)swzMalloc(sizeof(TMatrixHistogram)))) 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**)swzMalloc(cols*sizeof(int*)))) dw_Error(MEM_ERR);
for (j=cols-1; j >= 0; j--)
if (!(h->freq[i][j]=(int*)swzMalloc(intervals*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*)swzMalloc(cols*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*)swzMalloc(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--) swzFree(h->freq[i][j]);
swzFree(h->freq[i]);
}
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);
swzFree(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 *)swzMalloc(sizeof(TVectorHistogram))))
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*)swzMalloc(intervals*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);
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--) swzFree(h->freq[i]);
swzFree(h->freq);
swzFree(h->low);
swzFree(h->high);
FreeVector(h->Min);
FreeVector(h->Max);
swzFree(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 *)swzMalloc(sizeof(TScalarHistogram)))) dw_Error(MEM_ERR);
if (!(h->freq=(int*)swzMalloc(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)
{
swzFree(h->freq);
swzFree(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);
}

77
matlab/swz/c-code/utilities/DWCcode/histogram/dw_histogram.h Normal file
View File

@ -0,0 +1,77 @@
#ifndef __HISTOGRAMS__
#define __HISTOGRAMS__
#include "swzmatrix.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

10
mex/build/swz.am
View File

@ -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