simult_

 Simulates the model using a perturbation approach, given the path for the exogenous variables and the
 decision rules.

 INPUTS
    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number 
                        of auxilliary variables for lags and leads)
    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
    ex_      [double]   T*q matrix of innovations.
    iorder   [integer]  order of the taylor approximation.

 OUTPUTS
    y_       [double]   n*(T+1) time series for the endogenous variables.

 SPECIAL REQUIREMENTS
    none

0001 function y_=simult_(y0,dr,ex_,iorder)
0002 % Simulates the model using a perturbation approach, given the path for the exogenous variables and the
0003 % decision rules.
0004 %
0005 % INPUTS
0006 %    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number
0007 %                        of auxilliary variables for lags and leads)
0008 %    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
0009 %    ex_      [double]   T*q matrix of innovations.
0010 %    iorder   [integer]  order of the taylor approximation.
0011 %
0012 % OUTPUTS
0013 %    y_       [double]   n*(T+1) time series for the endogenous variables.
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %    none
0017 
0018 % Copyright (C) 2001-2011 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 global M_ options_
0036 
0037 iter = size(ex_,1);
0038 
0039 y_ = zeros(size(y0,1),iter+M_.maximum_lag);
0040 y_(:,1) = y0;
0041 
0042 % stoch_simul sets k_order_solver=1 if order=3, but does so only locally, so we
0043 % have to do it here also
0044 if options_.order == 3
0045     options_.k_order_solver = 1;
0046 end
0047 
0048 if ~options_.k_order_solver
0049     if iorder==1
0050         y_(:,1) = y_(:,1)-dr.ys;
0051     end
0052 end
0053 
0054 if options_.k_order_solver% Call dynare++ routines.
0055     ex_ = [zeros(1,M_.exo_nbr); ex_];
0056     switch options_.order
0057       case 1
0058         [err, y_] = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
0059                                   y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),...
0060                                   zeros(M_.endo_nbr,1),dr.g_1);
0061       case 2
0062         [err, y_] = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
0063                                   y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
0064                                   dr.g_1,dr.g_2);
0065       case 3
0066         [err, y_] = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ...
0067                                   y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed,dr.ys(dr.order_var),dr.g_0, ...
0068                                   dr.g_1,dr.g_2,dr.g_3);
0069       otherwise
0070         error(['order = ' int2str(order) ' isn''t supported'])
0071     end
0072     mexErrCheck('dynare_simul_', err);
0073     y_(dr.order_var,:) = y_;
0074 else
0075     if options_.block
0076         if M_.maximum_lag > 0
0077             k2 = dr.state_var;
0078         else
0079             k2 = [];
0080         end;
0081         order_var = 1:M_.endo_nbr;
0082         dr.order_var = order_var;
0083     else
0084         k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
0085         k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr;
0086         order_var = dr.order_var;
0087     end;
0088     
0089     switch iorder
0090       case 1
0091         if isempty(dr.ghu)% For (linearized) deterministic models.
0092             for i = 2:iter+M_.maximum_lag
0093                 yhat = y_(order_var(k2),i-1);
0094                 y_(order_var,i) = dr.ghx*yhat;
0095             end
0096         elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
0097             y_(dr.order_var,:) = dr.ghu*transpose(ex_);
0098         else
0099             epsilon = dr.ghu*transpose(ex_);
0100             for i = 2:iter+M_.maximum_lag
0101                 yhat = y_(order_var(k2),i-1);
0102                 y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
0103             end
0104         end
0105         y_ = bsxfun(@plus,y_,dr.ys);
0106       case 2
0107         constant = dr.ys(order_var)+.5*dr.ghs2;
0108         if options_.pruning
0109             y__ = y0;
0110             for i = 2:iter+M_.maximum_lag
0111                 yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
0112                 yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
0113                 epsilon = ex_(i-1,:)';
0114                 [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
0115                 mexErrCheck('A_times_B_kronecker_C', err);
0116                 [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
0117                 mexErrCheck('A_times_B_kronecker_C', err);
0118                 [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
0119                 mexErrCheck('A_times_B_kronecker_C', err);
0120                 y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
0121                     + abcOut1 + abcOut2 + abcOut3;
0122                 y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
0123             end
0124         else
0125             for i = 2:iter+M_.maximum_lag
0126                 yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
0127                 epsilon = ex_(i-1,:)';
0128                 [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
0129                 mexErrCheck('A_times_B_kronecker_C', err);
0130                 [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
0131                 mexErrCheck('A_times_B_kronecker_C', err);
0132                 [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
0133                 mexErrCheck('A_times_B_kronecker_C', err);
0134                 y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
0135                     + abcOut1 + abcOut2 + abcOut3;
0136             end
0137         end
0138     end
0139 end