simk

 function simk
 performs deterministic simulations with lead or lag on more than one
 period

 Currently used only for purely forward models.

 INPUTS
   ...
 OUTPUTS
   ...
 ALGORITHM
   Laffargue, Boucekkine, Juillard (LBJ)
   see Juillard (1996) Dynare: A program for the resolution and
   simulation of dynamic models with forward variables through the use
   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.

 SPECIAL REQUIREMENTS
   None.

0001 function simk
0002 % function simk
0003 % performs deterministic simulations with lead or lag on more than one
0004 % period
0005 %
0006 % Currently used only for purely forward models.
0007 %
0008 % INPUTS
0009 %   ...
0010 % OUTPUTS
0011 %   ...
0012 % ALGORITHM
0013 %   Laffargue, Boucekkine, Juillard (LBJ)
0014 %   see Juillard (1996) Dynare: A program for the resolution and
0015 %   simulation of dynamic models with forward variables through the use
0016 %   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.
0017 %
0018 % SPECIAL REQUIREMENTS
0019 %   None.
0020 %
0021 
0022 % Copyright (C) 1996-2012 Dynare Team
0023 %
0024 % This file is part of Dynare.
0025 %
0026 % Dynare is free software: you can redistribute it and/or modify
0027 % it under the terms of the GNU General Public License as published by
0028 % the Free Software Foundation, either version 3 of the License, or
0029 % (at your option) any later version.
0030 %
0031 % Dynare is distributed in the hope that it will be useful,
0032 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0034 % GNU General Public License for more details.
0035 %
0036 % You should have received a copy of the GNU General Public License
0037 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0038 
0039 global M_ options_ oo_
0040 
0041 nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
0042 ny = size(M_.lead_lag_incidence,2) ;
0043 icc1 = M_.lead_lag_incidence(nk,:) > 0;
0044 
0045 for i = 1:M_.maximum_lead -1
0046     icc1 = [M_.lead_lag_incidence(nk-i,:) | icc1(1,:); icc1] ;
0047 end
0048 
0049 icc1 = find(icc1') ;
0050 iy = M_.lead_lag_incidence > 0 ;
0051 isc = cumsum(sum(iy',1))' ;
0052 iyr0 = find(M_.lead_lag_incidence') ;
0053 ncc1 = size(icc1,1) ;
0054 ncc = ncc1 + 1 ;
0055 ncs = size(iyr0,1) ;
0056 
0057 ky = zeros(ny,nk) ;            % indices of variables at each lead or lag
0058 lky = zeros(nk,1) ;
0059 for i = 1:nk
0060     j = find(M_.lead_lag_incidence(i,:))' ;
0061     if isempty(j)
0062         lky(i) = 0;
0063     else
0064         lky(i) = size(j,1) ;
0065         ky(1:lky(i),i) = j ;
0066     end
0067 end
0068 
0069 jwc = find(iy(2:M_.maximum_endo_lead+1,:)') ; % indices of columns for
0070                                               % triangularization
0071                                               % as many rows as lags in model
0072 
0073 if isempty(jwc)
0074     jwc = 0 ;
0075     ljwc = 0 ;
0076     temp = icc1 ;
0077 else
0078     ljwc = size(jwc,1) ;          % length of each row in jwc
0079     temp = union(jwc,icc1,'rows') ;      % prepares next iteration
0080 end
0081 
0082 j1 = ky(1:lky(1),1) ;
0083 lj1 = lky(1) ;
0084 
0085 for i = 2:M_.maximum_endo_lag
0086     [j1,lj1] = ffill(j1,lj1,selif(temp+(i-1)*ny,temp <= ny)) ;
0087     if M_.maximum_endo_lead == 1
0088         if lky(i+M_.maximum_endo_lead) > 0
0089             [jwc,ljwc] = ffill(jwc,ljwc, ky(1:lky(i+M_.maximum_endo_lead),i+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny) ;
0090             if ljwc(i) == 0
0091                 temp = icc1;
0092             else
0093                 temp = union(jwc(1:ljwc(i),i),icc1,'rows') ;
0094             end
0095         else
0096             [jwc,ljwc] = ffill(jwc,ljwc,[]) ;
0097             temp = icc1 ;
0098         end
0099     else
0100         temp = temp(lj1(i)+1:size(temp,1),:) - ny ;
0101         if lky(i+M_.maximum_endo_lead) > 0
0102             [jwc,ljwc] = ffill(jwc,ljwc,[temp;ky(1:lky(i+M_.maximum_endo_lead),i+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny]);
0103         else
0104             [jwc,ljwc] = ffill(jwc,ljwc,temp) ;
0105         end
0106         temp = union(jwc(1:ljwc(i),i),icc1,'rows') ;
0107     end
0108 end
0109 
0110 [j1,lj1] = ffill(j1,lj1,selif(temp+M_.maximum_endo_lag*ny, temp <= ny)) ;
0111 ltemp = zeros(M_.maximum_endo_lag,1) ;
0112 jwc1 = zeros(ncc1,M_.maximum_endo_lag) ;
0113 
0114 for i = 1:M_.maximum_endo_lag
0115     temp = union(jwc(1:ljwc(i),i),icc1,'rows') ;
0116     ltemp(i) = size(temp,1) ;
0117     if ljwc(i) > 0
0118         jwc(1:ljwc(i),i) = indnv(jwc(1:ljwc(i),i),temp) ;
0119     end
0120     jwc1(:,i) = indnv(icc1,temp) ;
0121 end
0122 
0123 h1 = clock ;
0124 
0125 disp (['-----------------------------------------------------']) ;
0126 disp ('MODEL SIMULATION') ;
0127 fprintf ('\n') ;
0128 
0129 for iter = 1:options_.maxit_
0130     h2 = clock ;
0131     oo_.endo_simul = oo_.endo_simul(:);
0132     err_f = 0;
0133     
0134     fid = fopen([M_.fname '.swp'],'w+') ;
0135 
0136     it_ = 1+M_.maximum_lag ;
0137     ic = [1:ny] ;
0138     iyr = iyr0 ;
0139     i = M_.maximum_endo_lag+1 ;
0140     while (i>1) && (it_<=options_.periods+M_.maximum_endo_lag)
0141         h3 = clock ;
0142         [d1,jacobian] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params,oo_.steady_state, it_);
0143         d1 = -d1 ;
0144         err_f = max(err_f,max(abs(d1)));
0145         if lky(i) ~= 0
0146             j1i = ky(1:lky(i),i) ;
0147             w0 = jacobian(:,isc(i-1)+1:isc(i)) ;
0148         else
0149             w0 = [];
0150         end
0151         ttemp = iy(i+1:i+M_.maximum_endo_lead,:)' ;
0152         jwci = find(ttemp) ;
0153         if ~ isempty(jwci)
0154             w = jacobian(:,isc(i)+1:isc(i+M_.maximum_endo_lead)) ;
0155         end
0156         j = i ;
0157         while j <= M_.maximum_endo_lag
0158             if ~isempty(w0)
0159 
0160                 ofs = ((it_-M_.maximum_lag-M_.maximum_endo_lag+j-2)*ny)*ncc*8 ;
0161                 junk = fseek(fid,ofs,-1) ;
0162                 c = fread(fid,[ncc,ny],'float64')';
0163 
0164                 if isempty(jwci)
0165                     w = -w0*c(j1i,1:ncc1) ;
0166                     jwci = icc1 ;
0167                 else
0168                     iz = union(jwci,icc1,'rows') ;
0169                     ix = indnv(jwci,iz) ;
0170                     iy__ = indnv(icc1,iz) ;
0171                     temp = zeros(size(w,1),size(iz,1)) ;
0172                     temp(:,ix) = w;
0173                     temp(:,iy__) = temp(:,iy__)-w0*c(j1i,1:ncc1) ;
0174                     w = temp ;
0175                     jwci = iz ;
0176                     clear temp iz ix iy__ ;
0177                 end
0178                 d1 = d1-w0*c(j1i,ncc) ;
0179                 clear c ;
0180             end
0181             j = j + 1 ;
0182             if isempty(jwci)
0183                 j1i = [];
0184                 if lky(j+M_.maximum_endo_lead) ~= 0
0185                     jwci = ky(1:lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead) + (M_.maximum_endo_lead-1)*ny ;
0186                     w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
0187                 else
0188                     jwci = [] ;
0189                 end
0190             else
0191                 j1i = selif(jwci,jwci<(ny+1)) ;
0192                 w0 = w(:,1:size(j1i,1)) ;
0193                 if size(jwci,1) == size(j1i,1)
0194                     if lky(j+M_.maximum_endo_lead) ~= 0
0195                         jwci = ky(1:lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny ;
0196                         w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
0197                     else
0198                         jwci = [] ;
0199                     end
0200                 else
0201                     jwci = jwci(size(j1i,1)+1:size(jwci,1),:)-ny ;
0202                     w = w(:,size(j1i,1)+1:size(w,2)) ; 
0203                     if lky(j+M_.maximum_endo_lead) ~= 0
0204                         jwci = [ jwci; ky(1: lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny] ;
0205                         w = [w jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead))] ;
0206                         %         else
0207                         %           jwci = [] ;
0208                     end
0209                 end
0210             end
0211         end
0212         jwci = [indnv(jwci,icc1);ncc] ;
0213         w = [w d1] ;
0214         c = zeros(ny,ncc) ;
0215         c(:,jwci) = w0\w ;
0216         clear w w0 ;
0217 
0218         junk = fseek(fid,0,1) ;
0219         fwrite(fid,c','float64') ;
0220         clear c ;
0221 
0222         it_ = it_ + 1;
0223         ic = ic + ny ;
0224         iyr = iyr + ny ;
0225         i = i - 1 ;
0226     end
0227     icr0 = (it_-M_.maximum_lag-M_.maximum_endo_lag -1)*ny ;
0228     while it_ <= options_.periods+M_.maximum_lag
0229         [d1,jacobian] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params,oo_.steady_state, it_);
0230         d1 = -d1 ;
0231         err_f = max(err_f,max(abs(d1)));
0232         w0 = jacobian(:,1:isc(1)) ;
0233         w = jacobian(:,isc(1)+1:isc(1+M_.maximum_endo_lead)) ;
0234         j = 1 ;
0235         while j <= M_.maximum_endo_lag
0236             icr = j1(1:lj1(j),j)-(j-1)*ny ;
0237 
0238             ofs = ((icr0+(j-1)*ny+1)-1)*ncc*8 ;
0239             junk = fseek(fid,ofs,-1) ;
0240             c = fread(fid,[ncc,ny],'float64')' ;
0241 
0242             temp = zeros(ny,ltemp(j)) ;
0243             if ljwc(j) > 0
0244                 temp(:,jwc(1:ljwc(j),j)) = w ;
0245             end
0246             temp(:,jwc1(:,j))=temp(:,jwc1(:,j))-w0*c(icr,1:ncc1) ;
0247             w = temp ;
0248             clear temp ;
0249             d1 = d1-w0*c(icr,ncc) ;
0250             clear c ;
0251             j = j + 1 ;
0252             w0 = w(:,1:lj1(j)) ;
0253             if M_.maximum_endo_lead == 1
0254                 w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
0255             else
0256                 w = w(:,lj1(j)+1:size(w,2)) ;
0257 
0258                 if lky(j+M_.maximum_endo_lead) > 0
0259                     w = [w jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead))] ;
0260                 end
0261             end
0262         end
0263         c = w0\[w d1] ;
0264         d1 = [] ;
0265         clear w w0 ;
0266         junk = fseek(fid,0,1) ;
0267         fwrite(fid,c','float64') ;
0268         clear c ;
0269         it_ = it_ + 1 ;
0270         ic = ic + ny ;
0271         iyr = iyr + ny ;
0272         icr0 = icr0 + ny ;
0273     end
0274     if options_.terminal_condition == 1
0275 
0276         ofs = (((it_-M_.maximum_lag-2)*ny+1)-1)*ncc*8 ;
0277         junk = fseek(fid,ofs,-1) ;
0278         c = fread(fid,[ncc,ny],'float64')';
0279 
0280         for i = 1:M_.maximum_endo_lead
0281             w = tril(triu(ones(ny,ny+ncc1))) ;
0282             w(:,jwc1(:,M_.maximum_endo_lag)) = w(:,jwc1(:,M_.maximum_endo_lag))+c(:,1:ncc1) ;
0283             c = [w(:,ny+1:size(w,2))' c(:,ncc)]/w(:,1:ny) ;
0284 
0285             junk = fseek(fid,0,1) ;
0286             fwrite(fid,c','float64') ;
0287 
0288             it_ = it_+1 ;
0289             ic = ic + ny ;
0290 
0291         end
0292     end
0293     oo_.endo_simul = reshape(oo_.endo_simul,ny,options_.periods+M_.maximum_lag+M_.maximum_endo_lead) ;
0294     if options_.terminal_condition == 1
0295         hbacsup = clock ;
0296         c = bksupk(ny,fid,ncc,icc1) ;
0297         hbacsup = etime(clock,hbacsup) ;
0298         c = reshape(c,ny,options_.periods+M_.maximum_endo_lead)' ;
0299         y(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lead+M_.maximum_endo_lag)) = y(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lead+M_.maximum_endo_lag))+options_.slowc*c' ;
0300     else
0301         hbacsup = clock ;
0302         c = bksupk(ny,fid,ncc,icc1) ;
0303         hbacsup = etime(clock,hbacsup) ;
0304         c = reshape(c,ny,options_.periods)' ;
0305         oo_.endo_simul(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lag)) = oo_.endo_simul(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lag))+options_.slowc*c' ;
0306     end
0307 
0308     fclose(fid) ;
0309 
0310     h2 = etime(clock,h2) ;
0311     [junk,i1] = max(abs(c));
0312     [junk,i2] = max(junk);
0313     disp(['variable ' M_.endo_names(i2,:) ' period ' num2str(i1(i2))])
0314     err = max(max(abs(c./options_.scalv'))) ;
0315     disp ([num2str(iter) '-     err = ' num2str(err)]) ;
0316     disp (['err_f = ' num2str(err_f)])
0317     disp (['    Time of this iteration  : ' num2str(h2)]) ;
0318     if options_.timing
0319         disp (['        Back substitution               : ' num2str(hbacsup)]) ;
0320     end
0321     if err < options_.dynatol.f
0322         h1 = etime(clock,h1) ;
0323         fprintf ('\n') ;
0324         disp (['        Total time of simulation        : ' num2str(h1)]) ;
0325         fprintf ('\n') ;
0326         disp (['        Convergence achieved.']) ;
0327         disp (['-----------------------------------------------------']) ;
0328         fprintf ('\n') ;
0329         return ;
0330     end
0331 end
0332 disp(['WARNING : the maximum number of iterations is reached.']) ;
0333 fprintf ('\n') ;
0334 disp (['-----------------------------------------------------']) ;