diff --git a/matlab/resid.m b/matlab/resid.m
index f01fd5630..0c1c46990 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -75,7 +75,6 @@ function resid(period)
     end
     disp(['Residual for equation number ' int2str(i) ' is equal to ' num2str(z(i))])
   end  
-  
   for i = 1:2
     disp(' ')
   end
diff --git a/matlab/save_params_and_steady_state.m b/matlab/save_params_and_steady_state.m
index 6c3c29fe5..477ae0107 100644
--- a/matlab/save_params_and_steady_state.m
+++ b/matlab/save_params_and_steady_state.m
@@ -52,5 +52,4 @@ function save_params_and_steady_state(filename)
   for i = 1:M_.exo_nbr
     stored_values.(deblank(M_.exo_names(i,:))) = oo_.exo_steady_state(i);
   end
-  
-  save(filename, 'stored_values');
+  save('-v6',filename, 'stored_values');
diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc
index b50a06ca3..b9348f6d6 100644
--- a/preprocessor/BlockTriangular.cc
+++ b/preprocessor/BlockTriangular.cc
@@ -303,7 +303,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
-          //cout << "count_Block = " << count_Block << " i = " << i << " size_other_endo = " << tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] << "\n";
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
           ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
@@ -428,9 +427,10 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const
         free(ModelBlock->Block_List[blk].Variable);
         free(ModelBlock->Block_List[blk].Exogenous);
         free(ModelBlock->Block_List[blk].Own_Derivative);
+        free(ModelBlock->Block_List[blk].Other_Endogenous);
         for (i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
           {
-            if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size)
+            if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size*/)
               {
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
@@ -438,8 +438,13 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
+                free(ModelBlock->Block_List[blk].IM_lead_lag[i].u_other_endo);
+                free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_other_endo);
+                free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_other_endo);
+                free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index_other_endo);
+                free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index_other_endo);
               }
-            if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo)
+            if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo*/)
               {
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
                 free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
@@ -539,8 +544,6 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
           tmp_output.str("");
           lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
           tmp_s << "y(it_, " << Index_Var_IM[first_count_equ].index+1 << ")";
-          //cout << "tmp_s=" << tmp_s.str() << " tmp_output=" << tmp_output.str() << "  " << bool(tmp_output.str()==tmp_s.str()) << " " << BlockSim(Simulation_Type)
-          //     << " first_count_equ=" << first_count_equ << " equation=" << Index_Equ_IM[first_count_equ].index << "\n";
           //Determine whether the equation could be evaluated rather than to be solved
           if (tmp_output.str()==tmp_s.str())
             {
@@ -553,7 +556,6 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
             {
               tmp_output.str("");
               rhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms);
-              //cout << "sec tmp_s=" << tmp_s.str() << " tmp_output=" << tmp_output.str() << "  " << bool(tmp_output.str()==tmp_s.str()) << " " << BlockSim(Simulation_Type) << "\n";
               if (tmp_output.str()==tmp_s.str())
                 {
                   if (Simulation_Type==SOLVE_BACKWARD_SIMPLE)
@@ -568,10 +570,8 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
                    || ((prev_Type ==  EVALUATE_BACKWARD_R || prev_Type ==  EVALUATE_BACKWARD) && (Simulation_Type == EVALUATE_BACKWARD_R || Simulation_Type == EVALUATE_BACKWARD))
                  )
                 {
-                  //cout << "Type[0].first=" << Type[0].first << " Type[0].second= " << Type[0].second << "\n";
                   BlockSimulationType c_Type = (Type[Type.size()-1]).first;
                   int c_Size = (Type[Type.size()-1]).second;
-                  //cout << "i=" << i << " Type.size()=" << Type.size() << " c_Size=" << c_Size << "\n";
                   Type[Type.size()-1]=make_pair(c_Type, ++c_Size);
                 }
               else
@@ -582,11 +582,8 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
         }
       else
         {
-          /*for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
-            cout << Index_Equ_IM[count_equ ].index+1 <<  "  " << Index_Var_IM[count_equ ].index+1 << "\n";*/
           Type.push_back(make_pair(Simulation_Type, Blck_Size));
         }
-      //cout << "Type.size()= " << Type.size() << " BlockSim(Simulation_Type) = " << BlockSim(Simulation_Type) << "\n";
       prev_Type = Simulation_Type;
     }
   return(Type);
@@ -720,7 +717,6 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
   Nb_SimulBlocks = 0;
   for (t_type::const_iterator it = Type.begin(); it!=Type.end(); it++)
     {
-      //cout << "Block " << i++ << " Type=" << BlockSim(it->first) << " Size=" << it->second << "\n";
       if (it->first==SOLVE_FORWARD_COMPLETE || it->first==SOLVE_BACKWARD_COMPLETE || it->first==SOLVE_TWO_BOUNDARIES_COMPLETE)
         {
           Nb_SimulBlocks++;
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 51d907f5a..7f7d5e717 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -828,6 +828,37 @@ DynaTypeStatement::writeOutput(ostream &output, const string &basename) const
          << "',var_list_);" << endl;
 }
 
+SaveParamsAndSteadyStateStatement::SaveParamsAndSteadyStateStatement(const string &filename_arg) :
+  filename(filename_arg)
+{
+}
+
+void
+SaveParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
+{
+  output << "save_params_and_steady_state('" << filename << "');" << endl;
+}
+
+LoadParamsAndSteadyStateStatement::LoadParamsAndSteadyStateStatement(const string &filename_arg) :
+  filename(filename_arg)
+{
+}
+
+void
+LoadParamsAndSteadyStateStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.load_params_and_steady_state_present = true;
+  mod_file_struct.load_params_and_steady_state_filename = filename;
+}
+
+
+void
+LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
+{
+  output << "load_params_and_steady_state('" << filename << "');" << endl;
+}
+
+
 ModelComparisonStatement::ModelComparisonStatement(const filename_list_type &filename_list_arg,
                                                    const OptionsList &options_list_arg) :
   filename_list(filename_list_arg),
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index e3fb2760c..a294e3393 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -96,7 +96,7 @@ class ParsingDriver;
 %token <string_val> INT_NUMBER
 %token INV_GAMMA1_PDF INV_GAMMA2_PDF IRF
 %token KALMAN_ALGO KALMAN_TOL
-%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR LU
+%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LU
 %token MARKOWITZ MARGINAL_DENSITY MAX
 %token METHOD MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN
 %token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS
@@ -109,7 +109,7 @@ class ParsingDriver;
 %token PRINT PRIOR_TRUNC PRIOR_ANALYSIS POSTERIOR_ANALYSIS
 %token <string_val> QUOTED_STRING
 %token QZ_CRITERIUM
-%token RELATIVE_IRF REPLIC RPLOT
+%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
 %token SHOCKS SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER SOLVE_ALGO
 %token SPARSE SPARSE_DLL STDERR STEADY STOCH_SIMUL
 %token TEX RAMSEY_POLICY PLANNER_DISCOUNT
@@ -194,6 +194,8 @@ statement : declaration
           | dynare_sensitivity
           | homotopy_setup
           | forecast
+          | load_params_and_steady_state
+          | save_params_and_steady_state
           ;
 
 declaration : parameters
@@ -1088,6 +1090,16 @@ dynasave : DYNASAVE '(' filename ')'';'
            { driver.run_dynasave($3); }
          ;
 
+load_params_and_steady_state: LOAD_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
+          {driver.run_load_params_and_steady_state($3);}
+          ;
+
+save_params_and_steady_state: SAVE_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
+          {driver.run_save_params_and_steady_state($3);}
+          ;
+
+
+
 model_comparison : MODEL_COMPARISON mc_filename_list ';'
                    { driver.run_model_comparison(); }
                  | MODEL_COMPARISON '(' o_marginal_density ')' mc_filename_list ';'
@@ -1234,9 +1246,9 @@ dynare_sensitivity_option : o_gsa_identification
                           | o_mode_file
                           | o_gsa_trans_ident
                           ;
- 
 
-homotopy_setup: HOMOTOPY_SETUP ';' homotopy_list END 
+
+homotopy_setup: HOMOTOPY_SETUP ';' homotopy_list END
                { driver.end_homotopy();};
 
 homotopy_list : homotopy_item
@@ -1263,6 +1275,7 @@ forecast_option: o_periods
           | o_conf_sig
           ;
 
+
 number : INT_NUMBER
        | FLOAT_NUMBER
        ;
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 7def0fef8..980f2b670 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -29,7 +29,7 @@ using namespace std;
     DynareFlex::lex(Dynare::parser::semantic_type *yylval,     \
                     Dynare::parser::location_type *yylloc,     \
                     ParsingDriver &driver)
- 
+
 // Shortcut to access tokens defined by Bison
 typedef Dynare::parser::token token;
 
@@ -117,6 +117,9 @@ int sigma_e = 0;
 <INITIAL>dynasave 	{BEGIN DYNARE_STATEMENT; return token::DYNASAVE;}
 <INITIAL>model_comparison 	{BEGIN DYNARE_STATEMENT; return token::MODEL_COMPARISON;}
 
+<INITIAL>load_params_and_steady_state  {BEGIN DYNARE_STATEMENT; return token::LOAD_PARAMS_AND_STEADY_STATE;}
+<INITIAL>save_params_and_steady_state  {BEGIN DYNARE_STATEMENT; return token::SAVE_PARAMS_AND_STEADY_STATE;}
+
 <INITIAL>steady {BEGIN DYNARE_STATEMENT; return token::STEADY;}
 <INITIAL>check {BEGIN DYNARE_STATEMENT; return token::CHECK;}
 <INITIAL>simul {BEGIN DYNARE_STATEMENT; return token::SIMUL;}
diff --git a/preprocessor/Makefile b/preprocessor/Makefile
index 204754aa8..4dd2f09c1 100644
--- a/preprocessor/Makefile
+++ b/preprocessor/Makefile
@@ -2,7 +2,7 @@ include Makefile.include
 
 ifeq ($(shell uname -o), Cygwin)
 	DYNARE_M = dynare_m.exe
-else 
+else
 	DYNARE_M = dynare_m
 endif
 
@@ -25,6 +25,7 @@ OBJS = \
 	ParsingDriver.o \
 	DataTree.o \
 	ModFile.o \
+	MatlabFile.o \
 	Statement.o \
 	ExprNode.o \
 	ModelNormalization.o \
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index fa5cc04ad..554f80d76 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -37,13 +37,32 @@ ModFile::~ModFile()
 }
 
 
+
 void
 ModFile::evalAllExpressions()
 {
   //Evaluate Parameters
-  cout << "Evaluating expressions ...";
+
   InitParamStatement *it;
-  int j=0;
+  vector< vector<double> >::iterator it2;
+  ostringstream constant;
+  NodeID tmp_id;
+  CollectStruct collect_struct;
+  int j=0, k;
+  if(mod_file_struct.load_params_and_steady_state_present)
+    {
+      cout << "Reading " << mod_file_struct.load_params_and_steady_state_filename << " ...";
+      matlab_file.MatFileRead(mod_file_struct.load_params_and_steady_state_filename);
+      string sname="stored_values";
+      bool tmp_b=matlab_file.Collect(sname, collect_struct);
+      matlab_file.Delete();
+      if(!tmp_b)
+        {
+          cout << "The structure " << sname << " is not found in " << mod_file_struct.load_params_and_steady_state_filename << "\n";
+        }
+      cout << "done\n";
+    }
+  cout << "Evaluating expressions ...";
   for(vector<Statement *>::const_iterator it1=statements.begin();it1!=statements.end(); it1++)
     {
       it=dynamic_cast<InitParamStatement *>(*it1);
@@ -63,6 +82,23 @@ ModFile::evalAllExpressions()
            }
         }
     }
+  if(mod_file_struct.load_params_and_steady_state_present && j!=symbol_table.parameter_nbr)
+    {
+      //Reading a Mat-File
+      for(k=0;k <symbol_table.parameter_nbr; k++)
+        {
+          if(global_eval_context.find(make_pair(k, eParameter))==global_eval_context.end())
+            {
+              map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.find(symbol_table.getNameByID(eParameter, k));
+              if(it2!=collect_struct.variable_double_name.end())
+                {
+                  j++;
+                  vector<double>::iterator it=it2->second.begin();
+                  global_eval_context[make_pair(k, eParameter)]=*it;
+                }
+            }
+        }
+    }
   if (j!=symbol_table.parameter_nbr)
     {
       cout << "Warning: Uninitialized parameters: \n";
@@ -71,7 +107,6 @@ ModFile::evalAllExpressions()
           if(global_eval_context.find(make_pair(j, eParameter))==global_eval_context.end())
             cout << " " << symbol_table.getNameByID(eParameter, j) << "\n";
         }
-
     }
   //Evaluate variables
   for(InitOrEndValStatement::init_values_type::const_iterator it=init_values.begin(); it!=init_values.end(); it++)
@@ -90,6 +125,82 @@ ModFile::evalAllExpressions()
           cout << "error in evaluation of variable\n";
         }
     }
+  if(mod_file_struct.load_params_and_steady_state_present && init_values.size()<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
+    {
+      for(j=0;j <symbol_table.endo_nbr; j++)
+        {
+          if(global_eval_context.find(make_pair(j, eEndogenous))==global_eval_context.end())
+            {
+              //it2=mat_file.variable.find(symbol_table.getNameByID(eEndogenous, j));
+              map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.find(symbol_table.getNameByID(eEndogenous, j));
+              if(it2!=collect_struct.variable_double_name.end())
+                {
+                  vector<double>::iterator it=it2->second.begin();
+                  global_eval_context[make_pair(j, eEndogenous)]=*it;
+                  constant.str("");
+                  if(*it>=0)
+                    {
+                      constant << *it;
+                      tmp_id=expressions_tree.AddNumConstant(constant.str());
+                    }
+                  else
+                    {
+                      constant << -*it;
+                      tmp_id=expressions_tree.AddUMinus(expressions_tree.AddNumConstant(constant.str()));
+                    }
+                  init_values.push_back(make_pair(it2->first, tmp_id));
+                }
+            }
+        }
+      for(j=0;j <symbol_table.exo_nbr; j++)
+        {
+          if(global_eval_context.find(make_pair(j, eExogenous))==global_eval_context.end())
+            {
+              map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.find(symbol_table.getNameByID(eExogenous, j));
+              if(it2!=collect_struct.variable_double_name.end())
+                {
+                  vector<double>::iterator it=it2->second.begin();
+                  global_eval_context[make_pair(j, eExogenous)]=*it;
+                  constant.str("");
+                  if(*it>=0)
+                    {
+                      constant << *it;
+                      tmp_id=expressions_tree.AddNumConstant(constant.str());
+                    }
+                  else
+                    {
+                      constant << -*it;
+                      tmp_id=expressions_tree.AddUMinus(expressions_tree.AddNumConstant(constant.str()));
+                    }
+                  init_values.push_back(make_pair(it2->first, tmp_id));
+                }
+            }
+        }
+      for(j=0;j <symbol_table.exo_det_nbr; j++)
+        {
+          if(global_eval_context.find(make_pair(j, eExogenous))==global_eval_context.end())
+            {
+              map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.find(symbol_table.getNameByID(eExogenous, j));
+              if(it2!=collect_struct.variable_double_name.end())
+                {
+                  vector<double>::iterator it=it2->second.begin();
+                  global_eval_context[make_pair(j, eExogenous)]=*it;
+                  constant.str("");
+                  if(*it>=0)
+                    {
+                      constant << *it;
+                      tmp_id=expressions_tree.AddNumConstant(constant.str());
+                    }
+                  else
+                    {
+                      constant << -*it;
+                      tmp_id=expressions_tree.AddUMinus(expressions_tree.AddNumConstant(constant.str()));
+                    }
+                  init_values.push_back(make_pair(it2->first, tmp_id));
+                }
+            }
+        }
+    }
   if(init_values.size()<symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr)
     {
       cout << "\nWarning: Uninitialized variable: \n";
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index aa94d1287..127d70578 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -628,7 +628,6 @@ end:
                     << ", equation=" << eq+1 << endl;
                   }
               }
-            //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
             for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
               {
                 k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -677,7 +676,6 @@ end:
                     << ", equation=" << eq+1 << endl;
                   }
               }
-            /*jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr;*/
             for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
               {
                 k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -694,7 +692,6 @@ end:
                     << ", equation=" << eq+1 << endl;
                   }
               }
-            //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
             for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
               {
                 k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -726,7 +723,6 @@ end:
                 int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
                 int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
                 int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                //Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << ", " << varr+1 << ")*y(it_, " << var+1 << ")";
                 output << "    g1(" << eqr+1 << ", " << varr+1 << ") = ";
                 writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
                 output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
@@ -833,7 +829,6 @@ end:
                     << ", equation=" << eq+1 << endl;
                   }
               }
-            //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
             for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
               {
                 k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -854,7 +849,6 @@ end:
                       }
                   }
               }
-            //output << "  else" << endl;
             output << "      varargout{1}=g1_x;\n";
             output << "      varargout{2}=g1_o;\n";
             output << "    end;\n";
@@ -955,7 +949,6 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
         if (ModelBlock->Block_List[j].Temporary_terms->size())
           output << "  " << sps << "% //Temporary variables" << endl;
         i=0;
-        //temporary_terms_type tt2;
         for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin();
              it != ModelBlock->Block_List[j].Temporary_terms->end(); it++)
           {
@@ -971,18 +964,14 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
         // The equations
         for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
           {
-            //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
             string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ;
             output << sps << "  % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : "
             << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
-            /*if (!lhs_rhs_done)
-              {*/
-                eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
-                lhs = eq_node->arg1;
-                rhs = eq_node->arg2;
-                tmp_output.str("");
-                lhs->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms);
-              /*}*/
+            eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
+            lhs = eq_node->arg1;
+            rhs = eq_node->arg2;
+            tmp_output.str("");
+            lhs->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms);
             output << "  ";
             switch (ModelBlock->Block_List[j].Simulation_Type)
               {
@@ -1006,7 +995,6 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str
               case SOLVE_FORWARD_COMPLETE:
               case SOLVE_TWO_BOUNDARIES_COMPLETE:
               case SOLVE_TWO_BOUNDARIES_SIMPLE:
-                //Uf[ModelBlock->Block_List[j].Equation[i]] << "b(" << i+1 << ") =  residual(" << i+1 << ")";
                 goto end;
               default:
 end:
@@ -1054,12 +1042,6 @@ end:
                     int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
                     int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
                     int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                    /*if(!IM[eqr*ModelBlock->Block_List[j].Size+varr])
-                      {
-                        Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1
-                                                                    << ", " << varr+1 << ")*y( " << var+1 << ")";
-                        IM[eqr*ModelBlock->Block_List[j].Size+varr]=true;
-                      }*/
                     output << "  g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + ";
                     writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
                     output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
@@ -1270,11 +1252,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
                     lhs->compile(code_file,false, output_type, temporary_terms, map_idx);
                     rhs->compile(code_file,true, output_type, temporary_terms, map_idx);
                     break;
-                    /*case SOLVE_TWO_BOUNDARIES_SIMPLE:
-                      v=ModelBlock->Block_List[j].Equation[i];
-                      Uf[v].eqr=i;
-                      Uf[v].Ufl=NULL;
-                      goto end;*/
                   case SOLVE_BACKWARD_COMPLETE:
                   case SOLVE_FORWARD_COMPLETE:
                     v=ModelBlock->Block_List[j].Equation[i];
@@ -1951,15 +1928,9 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
     mStaticModelFile << "  y_kmax=M_.maximum_lead;\n";
     mStaticModelFile << "  y_size=M_.endo_nbr;\n";
 
-    /*tmp_output.str("");
-    writeModelLocalVariables(tmp_output, oMatlabDynamicModel);
-    if (tmp_output.str().length()>0)
-      mStaticModelFile << tmp_output.str() << "\n";*/
-
 
     mStaticModelFile << "  if(length(varargin)>0)\n";
     mStaticModelFile << "    %A simple evaluation of the static model\n";
-    //mStaticModelFile << "    global it_;\n";
     mStaticModelFile << "    y=varargin{1}(:);\n";
     mStaticModelFile << "    ys=y;\n";
     mStaticModelFile << "    g1=[];\n";
@@ -2762,47 +2733,6 @@ ModelTree::writeOutput(ostream &output) const
                 if (it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end())
                   exogenous.push_back(ii);
               }
-            /*if ((block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
-              ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
-              ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
-              ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
-              && j+Block_size<(block_triangular.ModelBlock->Size))
-              {
-                bool OK=true;
-                evaluate=true;
-                while(j+Block_size<(block_triangular.ModelBlock->Size) && OK)
-                  {
-                    if(BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j].Simulation_Type)!=BlockTriangular::BlockSim(block_triangular.ModelBlock->Block_List[j+Block_size].Simulation_Type))
-                      OK=false;
-                    else
-                      {
-                        if(max_lag <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag )
-                          max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ;
-                        if(max_lead<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead)
-                          max_lead=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead;
-                        if(max_lag_endo <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Endo )
-                          max_lag_endo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Endo ;
-                        if(max_lead_endo<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Endo)
-                          max_lead_endo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Endo;
-                        if(max_lag_exo <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo )
-                          max_lag_exo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo ;
-                        if(max_lead_exo<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Exo)
-                          max_lead_exo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Exo;
-                        for(int i=0;i<block_triangular.ModelBlock->Block_List[j+Block_size].Size;i++)
-                          {
-                            tmp_s << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Variable[i]+1;
-                            tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Equation[i]+1;
-                          }
-                        for(int i=0;i<block_triangular.ModelBlock->Block_List[j+Block_size].nb_exo;i++)
-                          {
-                            int ii=block_triangular.ModelBlock->Block_List[j+Block_size].Exogenous[i];
-                            if(it_exogenous==exogenous.end())
-                              exogenous.push_back(ii);
-                          }
-                        Block_size+=block_triangular.ModelBlock->Block_List[j+Block_size].Size;
-                      }
-                  }
-              }*/
             output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n";
             output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n";
             output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n";
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index d8652cec7..1c106e7be 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -1007,6 +1007,24 @@ ParsingDriver::run_dynasave(string *filename)
   delete filename;
 }
 
+
+void
+ParsingDriver::run_load_params_and_steady_state(string *filename)
+{
+  mod_file->addStatement(new LoadParamsAndSteadyStateStatement(*filename));
+
+  delete filename;
+}
+
+void
+ParsingDriver::run_save_params_and_steady_state(string *filename)
+{
+  mod_file->addStatement(new SaveParamsAndSteadyStateStatement(*filename));
+  delete filename;
+}
+
+
+
 void
 ParsingDriver::add_mc_filename(string *filename, string *prior)
 {
diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc
index 17adf591b..1032af8a9 100644
--- a/preprocessor/Statement.cc
+++ b/preprocessor/Statement.cc
@@ -29,7 +29,8 @@ ModFileStructure::ModFileStructure() :
   ramsey_policy_present(false),
   order_option(0),
   bvar_density_present(false),
-  bvar_forecast_present(false)
+  bvar_forecast_present(false),
+  load_params_and_steady_state_present(false)
 {
 }
 
diff --git a/preprocessor/include/ComputingTasks.hh b/preprocessor/include/ComputingTasks.hh
index c9ba68e4b..27a71994b 100644
--- a/preprocessor/include/ComputingTasks.hh
+++ b/preprocessor/include/ComputingTasks.hh
@@ -287,6 +287,27 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
+class SaveParamsAndSteadyStateStatement : public Statement
+{
+private:
+  const string filename;
+public:
+  SaveParamsAndSteadyStateStatement(const string &filename_arg);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
+class LoadParamsAndSteadyStateStatement : public Statement
+{
+private:
+  const string filename;
+public:
+  LoadParamsAndSteadyStateStatement(const string &filename_arg);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  string get_filename() const {return(filename);};
+};
+
+
 class ModelComparisonStatement : public Statement
 {
 public:
diff --git a/preprocessor/include/ModFile.hh b/preprocessor/include/ModFile.hh
index 3dfea1c96..3e0ed90c6 100644
--- a/preprocessor/include/ModFile.hh
+++ b/preprocessor/include/ModFile.hh
@@ -30,6 +30,7 @@ using namespace std;
 #include "ModelTree.hh"
 #include "VariableTable.hh"
 #include "Statement.hh"
+#include "MatlabFile.hh"
 
 //! The abstract representation of a "mod" file
 class ModFile
@@ -45,6 +46,8 @@ public:
   DataTree expressions_tree;
   //! Model equations and their derivatives
   ModelTree model_tree;
+  //! MatFile reading
+  MatlabFile matlab_file;
   //! Option linear
   bool linear;
   //! Global evaluation context
diff --git a/preprocessor/include/ParsingDriver.hh b/preprocessor/include/ParsingDriver.hh
index 67faa0700..38ea75a45 100644
--- a/preprocessor/include/ParsingDriver.hh
+++ b/preprocessor/include/ParsingDriver.hh
@@ -315,6 +315,8 @@ public:
   void run_calib(int covar);
   void run_dynasave(string *filename);
   void run_dynatype(string *filename);
+  void run_load_params_and_steady_state(string *filename);
+  void run_save_params_and_steady_state(string *filename);
   void add_mc_filename(string *filename, string *prior = new string("1"));
   void run_model_comparison();
   //! Begin a planner_objective statement
diff --git a/preprocessor/include/Statement.hh b/preprocessor/include/Statement.hh
index d86129f7a..324e0a143 100644
--- a/preprocessor/include/Statement.hh
+++ b/preprocessor/include/Statement.hh
@@ -54,6 +54,10 @@ public:
   bool bvar_density_present;
   //! Whether a bvar_forecast statement is present
   bool bvar_forecast_present;
+  //! Wether load_params_and_steady_state is present
+  bool load_params_and_steady_state_present;
+  //! save the load_params_and_steady state_filename
+  string load_params_and_steady_state_filename;
 };
 
 class Statement