extended_path_parfor

 Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
 series of size T  is obtained by solving T perfect foresight models.

 INPUTS
  o initial_conditions     [double]    m*nlags array, where m is the number of endogenous variables in the model and
                                       nlags is the maximum number of lags.
  o sample_size            [integer]   scalar, size of the sample to be simulated.

 OUTPUTS
  o time_series            [double]    m*sample_size array, the simulations.

 ALGORITHM

 SPECIAL REQUIREMENTS

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.max_lag = M_.maximum_endo_lag;
0042 pfm.nyp = nnz(pfm.lead_lag_incidence(1,:));
0043 pfm.iyp = find(pfm.lead_lag_incidence(1,:)>0);
0044 pfm.ny0 = nnz(pfm.lead_lag_incidence(2,:));
0045 pfm.iy0 = find(pfm.lead_lag_incidence(2,:)>0);
0046 pfm.nyf = nnz(pfm.lead_lag_incidence(3,:));
0047 pfm.iyf = find(pfm.lead_lag_incidence(3,:)>0);
0048 pfm.nd = pfm.nyp+pfm.ny0+pfm.nyf;
0049 pfm.nrc = pfm.nyf+1;
0050 pfm.isp = [1:pfm.nyp];
0051 pfm.is = [pfm.nyp+1:pfm.ny+pfm.nyp];
0052 pfm.isf = pfm.iyf+pfm.nyp;
0053 pfm.isf1 = [pfm.nyp+pfm.ny+1:pfm.nyf+pfm.nyp+pfm.ny+1];
0054 pfm.iz = [1:pfm.ny+pfm.nyp+pfm.nyf];
0055 pfm.periods = options_.ep.periods;
0056 pfm.steady_state = oo_.steady_state;
0057 pfm.params = M_.params;
0058 pfm.i_cols_1 = nonzeros(pfm.lead_lag_incidence(2:3,:)');
0059 pfm.i_cols_A1 = find(pfm.lead_lag_incidence(2:3,:)');
0060 pfm.i_cols_T = nonzeros(pfm.lead_lag_incidence(1:2,:)');
0061 pfm.i_cols_j = 1:pfm.nd;
0062 pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
0063 pfm.dynamic_model = str2func([M_.fname,'_dynamic']);
0064 pfm.verbose = options_.ep.verbosity;
0065 pfm.maxit_ = options_.maxit_;
0066 pfm.tolerance = options_.dynatol.f;
0067 
0068 exo_nbr = M_.exo_nbr;
0069 periods = options_.periods;
0070 ep = options_.ep;
0071 steady_state = oo_.steady_state;
0072 dynatol = options_.dynatol;
0073 
0074 % Set default initial conditions.
0075 if isempty(initial_conditions)
0076     initial_conditions = oo_.steady_state;
0077 end
0078 
0079 % Set maximum number of iterations for the deterministic solver.
0080 options_.maxit_ = options_.ep.maxit;
0081 
0082 % Set the number of periods for the perfect foresight model
0083 periods = options_.ep.periods;
0084 pfm.periods = options_.ep.periods;
0085 pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
0086 
0087 % Set the algorithm for the perfect foresight solver
0088 options_.stack_solve_algo = options_.ep.stack_solve_algo;
0089 
0090 % Set check_stability flag
0091 do_not_check_stability_flag = ~options_.ep.check_stability;
0092 
0093 
0094 % Compute the first order reduced form if needed.
0095 %
0096 % REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
0097 % all the globals in a mat file called linear_reduced_form.mat;
0098 
0099 dr = struct();
0100 if options_.ep.init
0101     options_.order = 1;
0102     [dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
0103 end
0104 
0105 % Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
0106 options_.minimal_solving_period = 100;%options_.ep.periods;
0107 
0108 % Initialize the exogenous variables.
0109 make_ex_;
0110 
0111 % Initialize the endogenous variables.
0112 make_y_;
0113 
0114 % Initialize the output array.
0115 time_series = zeros(M_.endo_nbr,sample_size);
0116 
0117 % Set the covariance matrix of the structural innovations.
0118 variances = diag(M_.Sigma_e);
0119 positive_var_indx = find(variances>0);
0120 effective_number_of_shocks = length(positive_var_indx);
0121 stdd = sqrt(variances(positive_var_indx));
0122 covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
0123 covariance_matrix_upper_cholesky = chol(covariance_matrix);
0124 
0125 % Set seed.
0126 if options_.ep.set_dynare_seed_to_default
0127     set_dynare_seed('default');
0128 end
0129 
0130 % Set bytecode flag
0131 bytecode_flag = options_.ep.use_bytecode;
0132 
0133 % Simulate shocks.
0134 switch options_.ep.innovation_distribution
0135   case 'gaussian'
0136       oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
0137   otherwise
0138     error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
0139 end
0140 
0141 % Set future shocks (Stochastic Extended Path approach)
0142 if options_.ep.stochastic.status
0143     switch options_.ep.stochastic.method
0144       case 'tensor'
0145         switch options_.ep.stochastic.ortpol
0146           case 'hermite'
0147             [r,w] = gauss_hermite_weights_and_nodes(options_.ep.stochastic.nodes);
0148           otherwise
0149             error('extended_path:: Unknown orthogonal polynomial option!')
0150         end
0151         if options_.ep.stochastic.order*M_.exo_nbr>1
0152             for i=1:options_.ep.stochastic.order*M_.exo_nbr
0153                 rr(i) = {r};
0154                 ww(i) = {w};
0155             end
0156             rrr = cartesian_product_of_sets(rr{:});
0157             www = cartesian_product_of_sets(ww{:});
0158         else
0159             rrr = r;
0160             www = w;
0161         end
0162         www = prod(www,2);
0163         number_of_nodes = length(www);
0164         relative_weights = www/max(www);
0165         switch options_.ep.stochastic.pruned.status
0166           case 1
0167             jdx = find(relative_weights>options_.ep.stochastic.pruned.relative);
0168             www = www(jdx);
0169             www = www/sum(www);
0170             rrr = rrr(jdx,:);
0171           case 2
0172             jdx = find(weights>options_.ep.stochastic.pruned.level);
0173             www = www(jdx);
0174             www = www/sum(www);
0175             rrr = rrr(jdx,:);
0176           otherwise
0177             % Nothing to be done!
0178         end
0179         nnn = length(www);
0180       otherwise
0181         error('extended_path:: Unknown stochastic_method option!')
0182     end
0183 else
0184     rrr = zeros(1,effective_number_of_shocks);
0185     www = 1;
0186     nnn = 1;
0187 end
0188 
0189 % Initializes some variables.
0190 t  = 0;
0191 
0192 % Set waitbar (graphic or text  mode)
0193 hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
0194 set(hh,'Name','EP simulations.');
0195 
0196 if options_.ep.memory
0197     mArray1 = zeros(M_.endo_nbr,100,nnn,sample_size);
0198     mArray2 = zeros(M_.exo_nbr,100,nnn,sample_size);
0199 end
0200 
0201 % Main loop.
0202 while (t<sample_size)
0203     if ~mod(t,10)
0204         dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
0205     end
0206     % Set period index.
0207     t = t+1;
0208     shocks = oo_.ep.shocks(t,:);
0209     % Put it in oo_.exo_simul (second line).
0210     oo_.exo_simul(2,positive_var_indx) = shocks;
0211     parfor s = 1:nnn
0212         periods1 = periods;
0213         exo_simul_1 = zeros(periods1+2,exo_nbr);
0214         pfm1 = pfm;
0215         info_convergence = 0;
0216         switch ep.stochastic.ortpol
0217           case 'hermite'
0218             for u=1:ep.stochastic.order
0219                 exo_simul_1(2+u,positive_var_indx) = rrr(s,(((u-1)*effective_number_of_shocks)+1):(u*effective_number_of_shocks))*covariance_matrix_upper_cholesky;
0220             end
0221           otherwise
0222             error('extended_path:: Unknown orthogonal polynomial option!')
0223         end
0224         if ep.stochastic.order && ep.stochastic.scramble
0225             exo_simul_1(2+ep.stochastic.order+1:2+ep.stochastic.order+ep.stochastic.scramble,positive_var_indx) = ...
0226                 randn(ep.stochastic.scramble,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
0227         end
0228         if ep.init% Compute first order solution (Perturbation)...
0229             ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
0230             ex(1:size(exo_simul_1,1),:) = exo_simul_1;
0231             exo_simul_1 = ex;
0232             initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
0233             endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
0234         else 
0235             endo_simul_1 = repmat(steady_state,1,periods1+2);
0236         end
0237         % Solve a perfect foresight model.
0238         increase_periods = 0;
0239         endo_simul = endo_simul_1;
0240         while 1
0241             if ~increase_periods
0242                 if bytecode_flag
0243                     [flag,tmp] = bytecode('dynamic');
0244                 else
0245                     flag = 1;
0246                 end
0247                 if flag
0248                     [flag,tmp] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
0249                 end
0250                 info_convergence = ~flag;
0251             end
0252             if verbosity
0253                 if info_convergence
0254                     if t<10
0255                         disp(['Time:    ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0256                     elseif t<100
0257                         disp(['Time:   ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0258                     elseif t<1000
0259                         disp(['Time:  ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0260                     else
0261                         disp(['Time: ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
0262                     end
0263                 else
0264                     if t<10
0265                         disp(['Time:    ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0266                     elseif t<100
0267                         disp(['Time:   ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0268                     elseif t<1000
0269                         disp(['Time:  ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0270                     else
0271                         disp(['Time: ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
0272                     end
0273                 end
0274             end
0275             if do_not_check_stability_flag
0276                 % Exit from the while loop.
0277                 endo_simul_1 = tmp;
0278                 break
0279             else
0280                 % Test if periods is big enough.
0281                 % Increase the number of periods.
0282                 periods1 = periods1 + ep.step;
0283                 pfm1.periods = periods1;
0284                 pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
0285                 % Increment the counter.
0286                 increase_periods = increase_periods + 1;
0287                 if verbosity
0288                     if t<10
0289                         disp(['Time:    ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0290                     elseif t<100
0291                         disp(['Time:   ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
0292                     elseif t<1000
0293                         disp(['Time:  ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0294                     else
0295                         disp(['Time: ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
0296                     end
0297                 end
0298                 if info_convergence
0299                     % If the previous call to the perfect foresight model solver exited
0300                     % announcing that the routine converged, adapt the size of endo_simul_1
0301                     % and exo_simul_1.
0302                     endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
0303                     exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,size(shocks,2)) ];
0304                     tmp_old = tmp;
0305                 else
0306                     % If the previous call to the perfect foresight model solver exited
0307                     % announcing that the routine did not converge, then tmp=1... Maybe
0308                     % should change that, because in some circonstances it may usefull
0309                     % to know where the routine did stop, even if convergence was not
0310                     % achieved.
0311                     endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
0312                     exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,size(shocks,2)) ];
0313                 end
0314                 % Solve the perfect foresight model with an increased number of periods.
0315                 if bytecode_flag
0316                     [flag,tmp] = bytecode('dynamic');
0317                 else
0318                     flag = 1;
0319                 end
0320                 if flag
0321                     [flag,tmp] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
0322                 end
0323                 info_convergence = ~flag;
0324                 if info_convergence
0325                     % If the solver achieved convergence, check that simulated paths did not
0326                     % change during the first periods.
0327                     % Compute the maximum deviation between old path and new path over the
0328                     % first periods
0329                     delta = max(max(abs(tmp(:,2:ep.fp)-tmp_old(:,2:ep.fp))));
0330                     if delta < dynatol.x
0331                         % If the maximum deviation is close enough to zero, reset the number
0332                         % of periods to ep.periods
0333                         periods1 = ep.periods;
0334                         pfm1.periods = periods1;
0335                         pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
0336                         % Cut exo_simul_1 and endo_simul_1 consistently with the resetted
0337                         % number of periods and exit from the while loop.
0338                         exo_simul_1 = exo_simul_1(1:(periods1+2),:);
0339                         endo_simul_1 = endo_simul_1(:,1:(periods1+2));
0340                         break
0341                     end
0342                 else
0343                     % The solver did not converge... Try to solve the model again with a bigger
0344                     % number of periods, except if the number of periods has been increased more
0345                     % than 10 times.
0346                     if increase_periods==10;
0347                         if verbosity
0348                             if t<10
0349                                 disp(['Time:    ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0350                             elseif t<100
0351                                 disp(['Time:   ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0352                             elseif t<1000
0353                                 disp(['Time:  ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0354                             else
0355                                 disp(['Time: ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
0356                             end
0357                         end
0358                         % Exit from the while loop.
0359                         break
0360                     end
0361                 end% if info_convergence
0362             end
0363         end% while
0364         if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
0365             [INFO,tmp] = homotopic_steps(.5,.01,pfm1);
0366             if (~isstruct(INFO) && isnan(INFO)) || ~info_convergence
0367                 [INFO,tmp] = homotopic_steps(0,.01,pfm1);
0368                 if ~info_convergence
0369                     disp('Homotopy:: No convergence of the perfect foresight model solver!')
0370                     error('I am not able to simulate this model!');
0371                 else
0372                     info_convergence = 1;
0373                     endo_simul_1 = tmp;
0374                     if verbosity && info_convergence
0375                         disp('Homotopy:: Convergence of the perfect foresight model solver!')
0376                     end
0377                 end
0378             else
0379                 info_convergence = 1;
0380                 endo_simul_1 = tmp;
0381                 if verbosity && info_convergence
0382                     disp('Homotopy:: Convergence of the perfect foresight model solver!')
0383                 end
0384             end
0385         end
0386         % Save results of the perfect foresight model solver.
0387         if ep.memory
0388             mArray1(:,:,s,t) = endo_simul_1(:,1:100);
0389             mArrat2(:,:,s,t) = transpose(exo_simul_1(1:100,:));
0390         end
0391         results(:,s) = www(s)*endo_simul_1(:,2);
0392     end
0393     time_series(:,t) = sum(results,2);
0394     oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
0395     oo_.endo_simul(:,1) = time_series(:,t);
0396     oo_.endo_simul(:,end) = oo_.steady_state;
0397 end% (while) loop over t
0398 
0399 dyn_waitbar_close(hh);
0400 
0401 oo_.endo_simul = oo_.steady_state;
0402 
0403 if ep.memory
0404     save([M_.fname '_memory'],'mArray1','mArray2','www');
0405 end