diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m
index e261de89c3657bb5674c94cb7fb4d7c9f23b911d..2bef6c5a6b480a9bf812a98edcd681f884c38f3b 100644
--- a/matlab/perfect-foresight-models/det_cond_forecast.m
+++ b/matlab/perfect-foresight-models/det_cond_forecast.m
@@ -144,7 +144,8 @@ else
end
%Compute the initial path using the the steady-state
% steady-state
- for jj = 2 : (options_.periods + 2)
+ %for jj = 2 : (options_.periods + 2)
+ for jj = 2 : (range.ndat + 2)
oo_.endo_simul(:, jj) = oo_.steady_state;
end
missings = isnan(oo_.endo_simul(:,1));
@@ -161,6 +162,7 @@ else
options_.dynatol.f = 1e-7;
[Info, endo, exo] = bytecode('extended_path', plan, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods);
options_.dynatol.f = save_options_dynatol_f;
+
if Info == 0
oo_.endo_simul = endo;
oo_.exo_simul = exo;
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 48c60019015c0ed230465cd402bc05c650b29025..4b2dccb9be93a5f7db4712c7cf52b0101567e782 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -862,13 +862,13 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
double *y_save = (double *) mxMalloc(size_of_direction);
double *x_save = (double *) mxMalloc((periods + y_kmax + y_kmin) * col_x *sizeof(double));*/
int max_periods = max(periods, nb_periods);
- size_t size_of_direction = y_size*(max_periods + y_kmax + y_kmin)*sizeof(double);
+ size_t size_of_direction = y_size*(/*max_periods*/ nb_periods + y_kmax + y_kmin)*sizeof(double);
double *y_save = (double *) mxMalloc(size_of_direction);
test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction);
- double *x_save = (double *) mxMalloc((max_periods + y_kmax + y_kmin) * col_x *sizeof(double));
- test_mxMalloc(x_save, __LINE__, __FILE__, __func__, (max_periods + y_kmax + y_kmin) * col_x *sizeof(double));
+ double *x_save = (double *) mxMalloc((/*max_periods*/nb_periods + y_kmax + y_kmin) * col_x *sizeof(double));
+ test_mxMalloc(x_save, __LINE__, __FILE__, __func__, (/*max_periods*/nb_periods + y_kmax + y_kmin) * col_x *sizeof(double));
vector_table_conditional_local_type vector_table_conditional_local;
vector_table_conditional_local.clear();
@@ -878,11 +878,20 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
{
x_save[j] = x[j];
x[j] = 0;
- }
+ }
+ for (int j = 0; j < col_x; j++)
+ for (int i = nb_row_x; i < nb_periods; i++)
+ x_save[j * (nb_periods + y_kmin + y_kmax) + i] = 0;//x[(j * nb_row_x) + 1];
for (int j = 0; j < col_x; j++)
x[y_kmin + j * nb_row_x] = x_save[1 + y_kmin + j * nb_row_x];
- for (int i = 0; i < (y_size*(periods + y_kmax + y_kmin)); i++)
- y_save[i] = y[i];
+ //for (int i = 0; i < (y_size*(periods + y_kmax + y_kmin)); i++)
+ for (int j = 0; j < y_size; j++)
+ for (int i = 0; i < nb_periods; i++)
+ if (i < periods)
+ y_save[j + (i + y_kmin) * y_size] = y[ j + (i + y_kmin) * y_size];
+ else
+ y_save[j + (i + y_kmin) * y_size] = y[ j + (periods + y_kmin) * y_size];
+
if (endo_name_length_l < 8)
endo_name_length_l = 8;
bool old_print_it = print_it;
@@ -913,8 +922,11 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
mexPrintf("|%s|",elastic(dates[t], date_length+2, false).c_str());
mexEvalString("drawnow;");
}
- for (vector<s_plan>::const_iterator it = sextended_path.begin(); it != sextended_path.end(); it++)
- x[y_kmin + (it->exo_num - 1) * (periods + y_kmax + y_kmin)] = it->value[t];
+ for (vector<s_plan>::const_iterator it = sextended_path.begin(); it != sextended_path.end(); it++)
+ {
+ x[y_kmin + (it->exo_num - 1) * /*(nb_periods + y_kmax + y_kmin)*/ nb_row_x] = it->value[t];
+ }
+
it_code = Init_Code;
vector_table_conditional_local.clear();
if (table_conditional_global.size())
@@ -930,9 +942,10 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
y[j ] = y[ j + (y_kmin) * y_size];
}
for (int j = 0; j < col_x; j++)
- {
- x_save[t + y_kmin + j * nb_row_x] = x[y_kmin + j * nb_row_x];
- x[y_kmin + j * nb_row_x] = x_save[t + 1 + y_kmin + j * nb_row_x];
+ {
+ x_save[t + y_kmin + j * nb_row_x] = x[y_kmin + j * nb_row_x];
+ if (t < nb_periods)
+ x[y_kmin + j * nb_row_x] = x_save[t + 1 + y_kmin + j * nb_row_x];
}
if (old_print_it)
@@ -941,31 +954,31 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
for (unsigned int i = 0; i < endo_name_length; i++)
if (P_endo_names[CHAR_LENGTH*(max_res_idx+i*y_size)] != ' ')
res << P_endo_names[CHAR_LENGTH*(max_res_idx+i*y_size)];
- res1 << std::scientific << max_res;
+ res1 << std::scientific << max_res;
mexPrintf("%s|%s| %4d | x |\n",elastic(res.str(),endo_name_length_l+2, true).c_str(), elastic(res1.str(), real_max_length+2, false).c_str(), iter);
mexPrintf(line.c_str());
mexEvalString("drawnow;");
}
}
print_it = old_print_it;
- for (int j = 0; j < y_size; j++)
+ /*for (int j = 0; j < y_size; j++)
{
for(int k = nb_periods; k < periods; k++)
y_save[j + (k + y_kmin) * y_size] = y[ j + ( k - (nb_periods-1) + y_kmin) * y_size];
- }
- for (int i = 0; i < (y_size*(periods + y_kmax + y_kmin)); i++)
+ }*/
+ for (int i = 0; i < (y_size*(nb_periods + y_kmax + y_kmin)); i++)
y[i] = y_save[i];
for (int j = 0; j < col_x* nb_row_x; j++)
- x[j] = x_save[j];
+ x[j] = x_save[j];
if (Init_Code->second)
- mxFree(Init_Code->second);
+ mxFree(Init_Code->second);
if (y_save)
- mxFree(y_save);
+ mxFree(y_save);
if (x_save)
- mxFree(x_save);
+ mxFree(x_save);
nb_blocks = Block_Count+1;
if (T && !global_temporary_terms)
- mxFree(T);
+ mxFree(T);
return true;
}
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 8cd42ca09b4277690d6ee59402e2dbd29dfec6db..e999dd025a8d73ff83dbf2114d43bfb038a903f7 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -1129,7 +1129,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
{
if ((t == 0) && (var < (periods+y_kmax)*Size) && (lag == 0) && (vector_table_conditional_local.size()))
{
- flip_exo = vector_table_conditional_local[var].var_exo;
+ flip_exo = vector_table_conditional_local[var].var_exo;
#ifdef DEBUG
local_index = eq;
#endif
@@ -3408,9 +3408,11 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
{
int eq = index_vara[i+Size*(y_kmin)];
int flip_exo = vector_table_conditional_local[i].var_exo;
- double yy = -(res[i] + x[y_kmin + flip_exo*/*row_x*/nb_row_x]);
- direction[eq] = 0;
- x[flip_exo*/*row_x*/nb_row_x + y_kmin] += slowc_l * yy;
+ double yy = -(res[i] + x[y_kmin + flip_exo*nb_row_x]);
+ direction[eq] = 0;
+ //mexPrintf("solving for init x[%d]=%f\n", flip_exo*nb_row_x + y_kmin, x[flip_exo*nb_row_x + y_kmin]);
+ x[flip_exo*nb_row_x + y_kmin] += slowc_l * yy;
+ //mexPrintf(" ==> x[flip_exo=%d (%d)]=%f yy=%f, slowc_l=%f, res[%d]=%f, y[5]=%f\n", flip_exo, flip_exo*nb_row_x + y_kmin, x[flip_exo*nb_row_x + y_kmin], yy, slowc_l, i, res[i], y[5]);
}
else
{
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index eee79bf454770bcad0aedd67f98e7d51f2f38669..c67bc09207dd697160c0ad3d221f05340d53cf59 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -1080,9 +1080,32 @@ main(int nrhs, const char *prhs[])
{
y[i] = double (yd[i]);
ya[i] = double (yd[i]);
- }
- size_t y_size = row_y;
+ }
+ size_t y_size = row_y;
size_t nb_row_x = row_x;
+
+ if (extended_path)
+ {
+ if (max_periods > periods)
+ {
+ double* y_save = (double*) mxMalloc(y_size * (max_periods + y_kmin + y_kmax) * sizeof(double));
+ for (int j = 0; j < y_size; j++)
+ for (int i = 0; i < periods; i++)
+ y_save[j + (i + y_kmin) * y_size] = y[ j + (i + y_kmin) * y_size];
+
+ y = (double*) mxRealloc(y, y_size * (max_periods + y_kmin + y_kmax) * sizeof(double));
+ x = (double*) mxRealloc(x, nb_row_x * (max_periods + y_kmin + y_kmax) * sizeof(double));
+ for (int j = 0; j < nb_row_x * (max_periods + y_kmin + y_kmax); j++)
+ x[j] = 0;
+ for (int j = 0; j < y_size; j++)
+ for (int i = 0; i < max_periods; i++)
+ if (i < periods)
+ y[j + (i + y_kmin) * y_size] = y_save[ j + (i + y_kmin) * y_size];
+ else
+ y[j + (i + y_kmin) * y_size] = y_save[ j + (periods + y_kmin) * y_size];
+ mxFree(y_save);
+ }
+ }
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal,
markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms, steady_state,
@@ -1096,8 +1119,9 @@ main(int nrhs, const char *prhs[])
int nb_blocks = 0;
double *pind;
bool no_error = true;
+
if (extended_path)
- {
+ {
try
{
interprete.extended_path(f, f, evaluate, block, nb_blocks, max_periods, sextended_path, sconditional_extended_path, dates, table_conditional_global);
@@ -1150,16 +1174,29 @@ main(int nrhs, const char *prhs[])
pind[i] = residual[i];
}
else
- {
- plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
+ {
+ int out_periods;
+ if (extended_path)
+ out_periods = max_periods + y_kmin + y_kmax;
+ else
+ out_periods = row_y;
+ plhs[1] = mxCreateDoubleMatrix(/*int(row_y)*/out_periods, int(col_y), mxREAL);
pind = mxGetPr(plhs[1]);
- for (i = 0; i < row_y*col_y; i++)
- pind[i] = y[i];
+ for (i = 0; i < /*row_y*/out_periods*col_y; i++)
+ {
+ pind[i] = y[i];
+ }
+
}
}
else
- {
- plhs[1] = mxCreateDoubleMatrix(int(row_y), int(col_y), mxREAL);
+ {
+ int out_periods;
+ if (extended_path)
+ out_periods = max_periods;
+ else
+ out_periods = col_y;
+ plhs[1] = mxCreateDoubleMatrix(int(row_y), /*int(col_y)*/out_periods, mxREAL);
pind = mxGetPr(plhs[1]);
if (evaluate)
{
@@ -1168,8 +1205,10 @@ main(int nrhs, const char *prhs[])
pind[i] = residual[i];
}
else
- for (i = 0; i < row_y*col_y; i++)
- pind[i] = y[i];
+ for (i = 0; i < row_y*out_periods /** col_y*/; i++)
+ {
+ pind[i] = y[i];
+ }
}
if (nlhs > 2)
{
@@ -1223,11 +1262,14 @@ main(int nrhs, const char *prhs[])
}
}
else
- {
+ {
plhs[2] = mxCreateDoubleMatrix(int(row_x), int(col_x), mxREAL);
pind = mxGetPr(plhs[2]);
- for (i = 0; i < row_x*col_x; i++)
- pind[i] = x[i];
+ for (i = 0; i < row_x*col_x; i++)
+ {
+ pind[i] = x[i];
+ }
+
}
if (nlhs > 3)
{
@@ -1252,15 +1294,15 @@ main(int nrhs, const char *prhs[])
}
#else
Free_global();
-#endif
+#endif
if (x)
- mxFree(x);
+ mxFree(x);
if (y)
- mxFree(y);
+ mxFree(y);
if (ya)
- mxFree(ya);
+ mxFree(ya);
if (direction)
- mxFree(direction);
+ mxFree(direction);
#ifdef _MSC_VER_
/*fFreeResult =*/ FreeLibrary(hinstLib);
#endif