Description of extended

0001 function time_series = extended_path(initial_conditions,sample_size)
0002 % Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
0003 % series of size T  is obtained by solving T perfect foresight models.
0004 %
0005 % INPUTS
0006 %  o initial_conditions     [double]    m*nlags array, where m is the number of endogenous variables in the model and
0007 %                                       nlags is the maximum number of lags.
0008 %  o sample_size            [integer]   scalar, size of the sample to be simulated.
0009 %
0010 % OUTPUTS
0011 %  o time_series            [double]    m*sample_size array, the simulations.
0012 %
0013 % ALGORITHM
0014 %
0015 % SPECIAL REQUIREMENTS
0016 
0017 % Copyright (C) 2009, 2010, 2011 Dynare Team
0018 %
0019 % This file is part of Dynare.
0020 %
0021 % Dynare is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % Dynare is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0033 global M_ options_ oo_
0034 
0035 options_.verbosity = options_.ep.verbosity;
0036 verbosity = options_.ep.verbosity+options_.ep.debug;
0037 
0038 % Prepare a structure needed by the matlab implementation of the perfect foresight model solver
0039 pfm.lead_lag_incidence = M_.lead_lag_incidence;
0040 pfm.ny = M_.endo_nbr;
0041 pfm.Sigma_e = M_.Sigma_e;
0042 max_lag = M_.maximum_endo_lag;
0043 pfm.max_lag = max_lag;
0044 if pfm.max_lag > 0
0045     pfm.nyp = nnz(pfm.lead_lag_incidence(1,:));
0046     pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0);
0047 else
0048     pfm.nyp = 0;
0049     pfm.iyp = [];
0050 end
0051 pfm.ny0 = nnz(pfm.lead_lag_incidence(max_lag+1,:));
0052 pfm.iy0 = find(pfm.lead_lag_incidence(max_lag+1,:)>0);
0053 if M_.maximum_endo_lead
0054     pfm.nyf = nnz(pfm.lead_lag_incidence(max_lag+2,:));
0055     pfm.iyf = find(pfm.lead_lag_incidence(max_lag+2,:)>0);
0056 else
0057     pfm.nyf = 0;
0058     pfm.iyf = [];
0059 end
0060 pfm.nd = pfm.nyp+pfm.ny0+pfm.nyf;
0061 pfm.nrc = pfm.nyf+1;
0062 pfm.isp = [1:pfm.nyp];
0063 pfm.is = [pfm.nyp+1:pfm.ny+pfm.nyp];
0064 pfm.isf = pfm.iyf+pfm.nyp;
0065 pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1];
0066 pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf];
0067 pfm.periods = options_.ep.periods;
0068 pfm.steady_state = oo_.steady_state;
0069 pfm.params = M_.params;
0070 if M_.maximum_endo_lead
0071     pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(max_lag+(1:2),:)');
0072     pfm.i_cols_A1 = find(pfm.lead_lag_incidence(max_lag+(1:2),:)');
0073 else
0074     pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(max_lag+1,:)');
0075     pfm.i_cols_A1 = find(pfm.lead_lag_incidence(max_lag+1,:)');
0076 end
0077 if max_lag > 0
0078     pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1:2,:)');
0079 else
0080     pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1,:)');
0081 end
0082 pfm.i_cols_j = 1:pfm.nd;
0083 pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
0084 pfm.dynamic_model = str2func([M_.fname,'_dynamic']);
0085 pfm.verbose = options_.ep.verbosity;
0086 pfm.maxit_ = options_.maxit_;
0087 pfm.tolerance = options_.dynatol.f;
0088 
0089 exo_nbr = M_.exo_nbr;
0090 periods = options_.periods;
0091 ep = options_.ep;
0092 steady_state = oo_.steady_state;
0093 dynatol = options_.dynatol;
0094 
0095 % Set default initial conditions.
0096 if isempty(initial_conditions)
0097     initial_conditions = oo_.steady_state;
0098 end
0099 
0100 % Set maximum number of iterations for the deterministic solver.
0101 options_.maxit_ = options_.ep.maxit;
0102 
0103 % Set the number of periods for the perfect foresight model
0104 periods = options_.ep.periods;
0105 pfm.periods = options_.ep.periods;
0106 pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
0107 
0108 % keep a copy of pfm.i_upd
0109 i_upd = pfm.i_upd;
0110 
0111 % Set the algorithm for the perfect foresight solver
0112 options_.stack_solve_algo = options_.ep.stack_solve_algo;
0113 
0114 % Set check_stability flag
0115 do_not_check_stability_flag = ~options_.ep.check_stability;
0116 
0117 % Compute the first order reduced form if needed.
0118 %
0119 % REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
0120 % all the globals in a mat file called linear_reduced_form.mat;
0121 
0122 dr = struct();
0123 if options_.ep.init
0124     options_.order = 1;
0125     [dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
0126 end
0127 
0128 % Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
0129 options_.minimal_solving_period = 100;%options_.ep.periods;
0130 
0131 % Initialize the exogenous variables.
0132 make_ex_;
0133 
0134 % Initialize the endogenous variables.
0135 make_y_;
0136 
0137 % Initialize the output array.
0138 time_series = zeros(M_.endo_nbr,sample_size);
0139 
0140 % Set the covariance matrix of the structural innovations.
0141 variances = diag(M_.Sigma_e);
0142 positive_var_indx = find(variances>0);
0143 effective_number_of_shocks = length(positive_var_indx);
0144 stdd = sqrt(variances(positive_var_indx));
0145 covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
0146 covariance_matrix_upper_cholesky = chol(covariance_matrix);
0147 
0148 % (re)Set exo_nbr
0149 %exo_nbr = effective_number_of_shocks;
0150 
0151 % Set seed.
0152 if options_.ep.set_dynare_seed_to_default
0153     set_dynare_seed('default');
0154 end
0155 
0156 % Set bytecode flag
0157 bytecode_flag = options_.ep.use_bytecode;
0158 
0159 % Simulate shocks.
0160 switch options_.ep.innovation_distribution
0161   case 'gaussian'
0162       oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
0163   otherwise
0164     error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
0165 end
0166 
0167 % Initializes some variables.
0168 t  = 0;
0169 
0170 % Set waitbar (graphic or text  mode)
0171 hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
0172 set(hh,'Name','EP simulations.');
0173 
0174 % Main loop.
0175 while (t<sample_size)
0176     if ~mod(t,10)
0177         dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
0178     end
0179     % Set period index.
0180     t = t+1;
0181     shocks = oo_.ep.shocks(t,:);
0182     % Put it in oo_.exo_simul (second line).
0183     oo_.exo_simul(2,positive_var_indx) = shocks;
0184     periods1 = periods;
0185     exo_simul_1 = zeros(periods1+2,exo_nbr);
0186     exo_simul_1(2,:) = oo_.exo_simul(2,:);
0187     pfm1 = pfm;
0188     info_convergence = 0;
0189     if ep.init% Compute first order solution (Perturbation)...
0190         ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
0191         ex(1:size(exo_simul_1,1),:) = exo_simul_1;
0192         exo_simul_1 = ex;
0193         initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
0194         endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
0195     else
0196         if t==1
0197             endo_simul_1 = repmat(steady_state,1,periods1+2);
0198         end
0199     end
0200     % Solve a perfect foresight model.
0201     increase_periods = 0;
0202     % Keep a copy of endo_simul_1
0203     endo_simul = endo_simul_1;
0204     while 1
0205         if ~increase_periods
0206             if bytecode_flag && ~options_.ep.stochastic.order
0207                 [flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
0208             else
0209                 flag = 1;
0210             end
0211             if flag
0212                 if options_.ep.stochastic.order == 0
0213                     [flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
0214                 else
0215                     [flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
0216                 end
0217             end
0218             info_convergence = ~flag;
0219         end
0220         if verbosity
0221             if info_convergence
0222                 if t<10
0223                     disp(['Time:    ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0224                 elseif t<100
0225                     disp(['Time:   ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0226                 elseif t<1000
0227                     disp(['Time:  ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0228                 else
0229                     disp(['Time: ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0230                 end
0231             else
0232                 if t<10
0233                     disp(['Time:    ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0234                 elseif t<100
0235                     disp(['Time:   ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0236                 elseif t<1000
0237                     disp(['Time:  ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0238                 else
0239                     disp(['Time: ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0240                 end
0241             end
0242         end
0243         if do_not_check_stability_flag
0244             % Exit from the while loop.
0245             endo_simul_1 = tmp;
0246             break
0247         else
0248             % Test if periods is big enough.
0249             % Increase the number of periods.
0250             periods1 = periods1 + ep.step;
0251             pfm1.periods = periods1;
0252             pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
0253             % Increment the counter.
0254             increase_periods = increase_periods + 1;
0255             if verbosity
0256                 if t<10
0257                     disp(['Time:    ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0258                 elseif t<100
0259                     disp(['Time:   ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
0260                 elseif t<1000
0261                     disp(['Time:  ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0262                 else
0263                     disp(['Time: ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0264                 end
0265             end
0266             if info_convergence
0267                 % If the previous call to the perfect foresight model solver exited
0268                 % announcing that the routine converged, adapt the size of endo_simul_1
0269                 % and exo_simul_1.
0270                 endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
0271                 exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
0272                 tmp_old = tmp;
0273             else
0274                 % If the previous call to the perfect foresight model solver exited
0275                 % announcing that the routine did not converge, then tmp=1... Maybe
0276                 % should change that, because in some circonstances it may usefull
0277                 % to know where the routine did stop, even if convergence was not
0278                 % achieved.
0279                 endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
0280                 exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
0281             end
0282             % Solve the perfect foresight model with an increased number of periods.
0283             if bytecode_flag && ~options_.ep.stochastic.order
0284                 [flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
0285             else
0286                 flag = 1;
0287             end
0288             if flag
0289                 if options_.ep.stochastic.order == 0
0290                     [flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
0291                 else
0292                     [flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
0293                 end
0294             end
0295             info_convergence = ~flag;
0296             if info_convergence
0297                 % If the solver achieved convergence, check that simulated paths did not
0298                 % change during the first periods.
0299                 % Compute the maximum deviation between old path and new path over the
0300                 % first periods
0301                 delta = max(max(abs(tmp(:,2)-tmp_old(:,2))));
0302                 if delta < dynatol.x
0303                     % If the maximum deviation is close enough to zero, reset the number
0304                     % of periods to ep.periods
0305                     periods1 = ep.periods;
0306                     pfm1.periods = periods1;
0307                     pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
0308                     % Cut exo_simul_1 and endo_simul_1 consistently with the resetted
0309                     % number of periods and exit from the while loop.
0310                     exo_simul_1 = exo_simul_1(1:(periods1+2),:);
0311                     endo_simul_1 = endo_simul_1(:,1:(periods1+2));
0312                     break
0313                 end
0314             else
0315                 % The solver did not converge... Try to solve the model again with a bigger
0316                 % number of periods, except if the number of periods has been increased more
0317                 % than 10 times.
0318                 if increase_periods==10;
0319                     if verbosity
0320                         if t<10
0321                             disp(['Time:    ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0322                         elseif t<100
0323                             disp(['Time:   ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0324                         elseif t<1000
0325                             disp(['Time:  ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0326                         else
0327                             disp(['Time: ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0328                         end
0329                     end
0330                     % Exit from the while loop.
0331                     break
0332                 end
0333             end% if info_convergence
0334         end
0335     end% while
0336     if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
0337         if ~do_not_check_stability_flag
0338             periods1 = ep.periods;
0339             pfm1.periods = periods1;
0340             pfm1.i_upd = i_upd;
0341             exo_simul_1 = exo_simul_1(1:(periods1+2),:);
0342             endo_simul_1 = endo_simul_1(:,1:(periods1+2));
0343         end
0344         [INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,.5,.01,pfm1);
0345         if isstruct(INFO)
0346             info_convergence = INFO.convergence;
0347         else
0348             info_convergence = 0;
0349         end
0350         if ~info_convergence
0351             [INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,0,.01,pfm1);
0352             if isstruct(INFO)
0353                 info_convergence = INFO.convergence;
0354             else
0355                 info_convergence = 0;
0356             end
0357             if ~info_convergence
0358                 disp('Homotopy:: No convergence of the perfect foresight model solver!')
0359                 error('I am not able to simulate this model!');
0360             else
0361                 endo_simul_1 = tmp;
0362                 if verbosity && info_convergence
0363                     disp('Homotopy:: Convergence of the perfect foresight model solver!')
0364                 end
0365             end
0366         else
0367             info_convergence = 1;
0368             endo_simul_1 = tmp;
0369             if verbosity && info_convergence
0370                 disp('Homotopy:: Convergence of the perfect foresight model solver!')
0371             end
0372         end
0373     end
0374     % Save results of the perfect foresight model solver.
0375     time_series(:,t) = endo_simul_1(:,2);
0376     endo_simul_1(:,1:end-1) = endo_simul_1(:,2:end);
0377     endo_simul_1(:,1) = time_series(:,t);
0378     endo_simul_1(:,end) = oo_.steady_state;
0379 end% (while) loop over t
0380 
0381 dyn_waitbar_close(hh);
0382 
0383 oo_.endo_simul = oo_.steady_state;
extended_path

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
