diff --git a/matlab/disp_dr_sparse.m b/matlab/disp_dr_sparse.m new file mode 100644 index 000000000..bf5ce3f0f --- /dev/null +++ b/matlab/disp_dr_sparse.m @@ -0,0 +1,234 @@ +function disp_dr_sparse(dr,order,var_list) + +% Copyright (C) 2001 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + + global M_ + nx = 0; + nu = 0; + k = []; + klag = []; + k1 = []; + nspred = 0; + for i=1:length(M_.block_structure.block) + nspred = nspred + M_.block_structure.block(i).dr.nspred; + end; + ghu = zeros(M_.endo_nbr, M_.exo_nbr*(M_.maximum_exo_lag+M_.maximum_exo_lead+1)); + ghx = zeros(M_.endo_nbr, nspred); + for i=1:length(M_.block_structure.block) + nx = nx + size(M_.block_structure.block(i).dr.ghx,2); +% M_.block_structure.block(i).dr.ghx +% M_.block_structure.block(i).equation +% M_.block_structure.block(i).variable + ghx(M_.block_structure.block(i).equation, M_.block_structure.block(i).variable(find(M_.block_structure.block(i).lead_lag_incidence(1: M_.block_structure.block(i).maximum_endo_lag,:))) ) = M_.block_structure.block(i).dr.ghx; + if(M_.block_structure.block(i).exo_nbr) + nu = nu + size(M_.block_structure.block(i).dr.ghu,2); + ghu(M_.block_structure.block(i).equation, M_.block_structure.block(i).exogenous) = M_.block_structure.block(i).dr.ghu; + end + k_tmp = find(M_.block_structure.block(i).dr.kstate(:,2) <= M_.block_structure.block(i).maximum_lag+1); + k = [k ; k_tmp]; + klag = [klag ; M_.block_structure.block(i).dr.kstate(k_tmp,[1 2])]; + k1 = [k1 ; M_.block_structure.block(i).variable(M_.block_structure.block(i).dr.order_var)']; + end + nvar = size(var_list,1); + if nvar == 0 + nvar = length(k1); + ivar = [1:nvar]; + else + ivar=zeros(nvar,1); + for i=1:nvar + i_tmp = strmatch(var_list(i,:),M_.endo_names(k1,:),'exact'); + if isempty(i_tmp) + disp(var_list(i,:)); + error (['One of the variable specified does not exist']) ; + else + ivar(i) = i_tmp; + end + end + end + disp('POLICY AND TRANSITION FUNCTIONS') + % variable names + str = ' '; + for i=1:nvar + str = [str sprintf('%16s',M_.endo_names(k1(ivar(i)),:))]; + end + disp(str); + % + % constant + % + str = 'Constant '; + flag = 0; + for i=1:nvar + x = dr.ys(k1(ivar(i))); + if order > 1 + x = x + dr.ghs2(ivar(i))/2; + end + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + if order > 1 + str = '(correction) '; + flag = 0; + for i=1:nvar + x = dr.ghs2(ivar(i))/2; + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + % + % ghx + % + for k=1:nx + flag = 0; + str1 = sprintf('%s(%d)',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2); + str = sprintf('%-20s',str1); + for i=1:nvar + x = ghx(ivar(i),k); + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + % + % ghu + % + for k=1:nu + flag = 0; + str = sprintf('%-20s',M_.exo_names(k,:)); + for i=1:nvar + x = ghu(ivar(i),k); + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + + if order > 1 + % ghxx + for k = 1:nx + for j = 1:k + flag = 0; + str1 = sprintf('%s(%d),%s(%d)',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2, ... + M_.endo_names(k1(klag(j,1)),:),klag(j,2)-M_.maximum_lag-2); + str = sprintf('%-20s',str1); + for i=1:nvar + if k == j + x = dr.ghxx(ivar(i),(k-1)*nx+j)/2; + else + x = dr.ghxx(ivar(i),(k-1)*nx+j); + end + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + end + % + % ghuu + % + for k = 1:nu + for j = 1:k + flag = 0; + str = sprintf('%-20s',[M_.exo_names(k,:) ',' M_.exo_names(j,:)] ); + for i=1:nvar + if k == j + x = dr.ghuu(ivar(i),(k-1)*nu+j)/2; + else + x = dr.ghuu(ivar(i),(k-1)*nu+j); + end + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + end + % + % ghxu + % + for k = 1:nx + for j = 1:nu + flag = 0; + str1 = sprintf('%s(%d),%s',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2, ... + M_.exo_names(j,:)); + str = sprintf('%-20s',str1); + for i=1:nvar + x = dr.ghxu(ivar(i),(k-1)*nu+j); + if abs(x) > 1e-6 + flag = 1; + str = [str sprintf('%16.6f',x)]; + else + str = [str ' 0']; + end + end + if flag + disp(str) + end + end + end + end + +% $$$ dr.ghx +% $$$ dr.ghu +% $$$ dr.ghxx +% $$$ dr.ghuu +% $$$ dr.ghxu + +% 01/08/2001 MJ added test for order in printing quadratic terms +% 02/21/2001 MJ pass all variable names through deblank() +% 02/21/2001 MJ changed from f to g format to write numbers +% 10/09/2002 MJ corrected error on constant whith subset of variables + + + diff --git a/matlab/dr11_sparse.m b/matlab/dr11_sparse.m index 1a2f7cc05..2d9a981ef 100644 --- a/matlab/dr11_sparse.m +++ b/matlab/dr11_sparse.m @@ -17,7 +17,7 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - + %task info = 0; klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1; kstate = dr.kstate; @@ -46,7 +46,12 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi m = m+length(k); end if M_.exo_nbr & task~=1 + jacobia_ + jacobia_(:,nz+1:end) + b dr.ghu = -b\jacobia_(:,nz+1:end); + disp(['nz=' int2str(nz) ]); + dr.ghu end dr.eigval = eig(transition_matrix(dr,M_)); dr.rank = 0; @@ -156,12 +161,15 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi k1 = find(kstate(n4:nd,2) == M_.maximum_endo_lag+1); k2 = find(kstate(1:n3,2) == M_.maximum_endo_lag+2); + hx(k1,:) + gx(k2(nboth+1:end),:) dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)]; - + dr.ghx %lead variables actually present in the model j3 = nonzeros(kstate(:,3)); j4 = find(kstate(:,3)); % derivatives with respect to exogenous variables + disp(['M_.exo_nbr=' int2str(M_.exo_nbr)]); if M_.exo_nbr fu = aa(:,nz+(1:M_.exo_nbr)); a1 = b; @@ -216,7 +224,6 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:); end end - disp('end0'); if options_.order == 1 return end @@ -479,5 +486,4 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi end end - disp('end'); end \ No newline at end of file diff --git a/matlab/dr1_sparse.m b/matlab/dr1_sparse.m index acef2657e..43c39f3cc 100644 --- a/matlab/dr1_sparse.m +++ b/matlab/dr1_sparse.m @@ -210,23 +210,47 @@ function [dr,info,M_,options_,oo_] = dr1_sparse(dr,task,M_,options_,oo_) if options_.order == 1 [junk,jacobia_] = feval([M_.fname '_dynamic'],ones(M_.maximum_lag+M_.maximum_lead+1,1)*dr.ys',[oo_.exo_simul ... oo_.exo_det_simul], it_); + %full(jacobia_) dr.eigval = []; dr.nyf = 0; dr.rank = 0; + first_col_exo = M_.endo_nbr * (M_.maximum_endo_lag + M_.maximum_endo_lead + 1); for i=1:length(M_.block_structure.block) + %disp(['block = ' int2str(i)]); M_.block_structure.block(i).dr.Null=0; M_.block_structure.block(i).dr=set_state_space(M_.block_structure.block(i).dr,M_.block_structure.block(i)); col_selector=repmat(M_.block_structure.block(i).variable,1,M_.block_structure.block(i).maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead+1)+kron([M_.maximum_endo_lag-M_.block_structure.block(i).maximum_endo_lag:M_.maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead],M_.endo_nbr*ones(1,M_.block_structure.block(i).endo_nbr)); row_selector = M_.block_structure.block(i).equation; + %col_selector jcb_=jacobia_(row_selector,col_selector); - jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ; + jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ; + if M_.block_structure.block(i).exo_nbr>0 + col_selector = [ first_col_exo + ... + repmat(M_.block_structure.block(i).exogenous,1,M_.block_structure.block(i).maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead+1)+kron([M_.maximum_exo_lag-M_.block_structure.block(i).maximum_exo_lag:M_.maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead],M_.exo_nbr*ones(1,M_.block_structure.block(i).exo_nbr))]; + end + %col_selector + jcb_ = [ jcb_ jacobia_(row_selector,col_selector)]; + %full(jcb_) + hss_=0; %hessian(M_.block_structure.block(i).equation,M_.block_structure.block(i).variable); dra = M_.block_structure.block(i).dr; - M_.block_structure.block(i).exo_nbr=M_.exo_nbr; + %M_.block_structure.block(i).exo_nbr=M_.exo_nbr; [dra ,info,M_.block_structure.block(i),options_,oo_] = dr11_sparse(dra ,task,M_.block_structure.block(i),options_,oo_, jcb_, hss_); M_.block_structure.block(i).dr = dra; dr.eigval = [dr.eigval; dra.eigval]; - dr.nyf = dr.nyf + nnz(dra.kstate(:,2)>M_.block_structure.block(i).maximum_endo_lag+1); + nyf = nnz(dra.kstate(:,2)>M_.block_structure.block(i).maximum_endo_lag+1); + n_explod = nnz(abs(dra.eigval) > options_.qz_criterium); + if nyf ~= n_explod + disp(['EIGENVALUES in block ' int2str(i) ':']); + [m_lambda,ii]=sort(abs(dra.eigval)); + disp(sprintf('%16s %16s %16s\n','Modulus','Real','Imaginary')) + z=[m_lambda real(dra.eigval(ii)) imag(dra.eigval(ii))]'; + disp(sprintf('%16.4g %16.4g %16.4g\n',z)) + disp(['The rank condition is not satisfy in block ' int2str(i) ' :']); + disp([' ' int2str(nyf) ' forward-looking variable(s) for ' ... + int2str(n_explod) ' eigenvalue(s) larger than 1 in modulus']); + end + dr.nyf = dr.nyf + nyf; dr.rank = dr.rank + dra.rank; end; end diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe index f17104b74..bdcfdffa7 100755 Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ diff --git a/matlab/stoch_simul_sparse.m b/matlab/stoch_simul_sparse.m new file mode 100644 index 000000000..09efa3076 --- /dev/null +++ b/matlab/stoch_simul_sparse.m @@ -0,0 +1,295 @@ +function info=stoch_simul(var_list) + +% Copyright (C) 2001-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + + global M_ options_ oo_ it_ + + options_old = options_; + if options_.linear + options_.order = 1; + end + if (options_.order == 1) + options_.replic = 1; + end + + + TeX = options_.TeX; + + iter_ = max(options_.periods,1); + if M_.exo_nbr > 0 + oo_.exo_simul= ones(iter_ + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state'; + end + + check_model; + + [oo_.dr, info] = resol(oo_.steady_state,0); + + if info(1) + options_ = options_old; + print_info(info); + return + end + + oo_dr_kstate = []; + oo_dr_nstatic = 0; + for i=1:length(M_.block_structure.block) + oo_dr_kstate = [ oo_dr_kstate ; M_.block_structure.block(i).dr.kstate]; + oo_dr_nstatic = oo_dr_nstatic + M_.block_structure.block(i).dr.nstatic; + end + + if ~options_.noprint + disp(' ') + disp('MODEL SUMMARY') + disp(' ') + disp([' Number of variables: ' int2str(M_.endo_nbr)]) + disp([' Number of stochastic shocks: ' int2str(M_.exo_nbr)]) + disp([' Number of state variables: ' ... + int2str(length(find(oo_dr_kstate(:,2) <= M_.maximum_lag+1)))]) + disp([' Number of jumpers: ' ... + int2str(length(find(oo_dr_kstate(:,2) == M_.maximum_lag+2)))]) + disp([' Number of static variables: ' int2str(oo_dr_nstatic)]) + my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS'; + labels = deblank(M_.exo_names); + headers = strvcat('Variables',labels); + lh = size(labels,2)+2; + table(my_title,headers,labels,M_.Sigma_e,lh,10,6); + disp(' ') + disp_dr_sparse(oo_.dr,options_.order,var_list); + end + + if options_.simul == 0 & options_.nomoments == 0 + disp_th_moments(oo_.dr,var_list); + elseif options_.simul == 1 + if options_.periods == 0 + error('STOCH_SIMUL error: number of periods for the simulation isn''t specified') + end + if options_.periods < options_.drop + disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ... + ' than the number of observations to be DROPed']) + options_ =options_old; + return + end + oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr); + dyn2vec; + if options_.nomoments == 0 + disp_moments(oo_.endo_simul,var_list); + end + end + + + + if options_.irf + n = size(var_list,1); + if n == 0 + n = M_.endo_nbr; + ivar = [1:n]'; + var_list = M_.endo_names; + if TeX + var_listTeX = M_.endo_names_tex; + end + else + ivar=zeros(n,1); + if TeX + var_listTeX = []; + end + for i=1:n + i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact'); + if isempty(i_tmp) + error (['One of the specified variables does not exist']) ; + else + ivar(i) = i_tmp; + if TeX + var_listTeX = strvcat(var_listTeX,deblank(M_.endo_names_tex(i_tmp,:))); + end + end + end + end + if TeX + fidTeX = fopen([M_.fname '_IRF.TeX'],'w'); + fprintf(fidTeX,'%% TeX eps-loader file generated by stoch_simul.m (Dynare).\n'); + fprintf(fidTeX,['%% ' datestr(now,0) '\n']); + fprintf(fidTeX,' \n'); + end + olditer = iter_;% Est-ce vraiment utile ? Il y a la même ligne dans irf... + SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr); + cs = transpose(chol(SS)); + tit(M_.exo_names_orig_ord,:) = M_.exo_names; + if TeX + titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex; + end + for i=1:M_.exo_nbr + if SS(i,i) > 1e-13 + y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ... + options_.replic, options_.order); + if options_.relative_irf + y = 100*y/cs(i,i); + end + irfs = []; + mylist = []; + if TeX + mylistTeX = []; + end + for j = 1:n + assignin('base',[deblank(M_.endo_names(ivar(j),:)) '_' deblank(M_.exo_names(i,:))],... + y(ivar(j),:)'); + eval(['oo_.irfs.' deblank(M_.endo_names(ivar(j),:)) '_' ... + deblank(M_.exo_names(i,:)) ' = y(ivar(j),:);']); + if max(y(ivar(j),:)) - min(y(ivar(j),:)) > 1e-10 + irfs = cat(1,irfs,y(ivar(j),:)); + mylist = strvcat(mylist,deblank(var_list(j,:))); + if TeX + mylistTeX = strvcat(mylistTeX,deblank(var_listTeX(j,:))); + end + end + end + if options_.nograph == 0 + number_of_plots_to_draw = size(irfs,1); + [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw); + if nbplt == 0 + elseif nbplt == 1 + if options_.relative_irf + hh = figure('Name',['Relative response to' ... + ' orthogonalized shock to ' tit(i,:)]); + else + hh = figure('Name',['Orthogonalized shock to' ... + ' ' tit(i,:)]); + end + for j = 1:number_of_plots_to_draw + subplot(nr,nc,j); + plot(1:options_.irf,transpose(irfs(j,:)),'-k','linewidth',1); + hold on + plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + hold off + xlim([1 options_.irf]); + title(deblank(mylist(j,:)),'Interpreter','none'); + end + eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:))]); + saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) '.fig']); + end + if TeX + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for j = 1:number_of_plots_to_draw + fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist(j,:)),deblank(mylistTeX(j,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s}\n',M_.fname,deblank(tit(i,:))); + fprintf(fidTeX,'\\caption{Impulse response functions (orthogonalized shock to $%s$).}',titTeX(i,:)); + fprintf(fidTeX,'\\label{Fig:IRF:%s}\n',deblank(tit(i,:))); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end + % close(hh) + else + for fig = 1:nbplt-1 + if options_.relative_irf == 1 + hh = figure('Name',['Relative response to orthogonalized shock' ... + ' to ' tit(i,:) ' figure ' int2str(fig)]); + else + hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ... + ' figure ' int2str(fig)]); + end + for plt = 1:nstar + subplot(nr,nc,plt); + plot(1:options_.irf,transpose(irfs((fig-1)*nstar+plt,:)),'-k','linewidth',1); + hold on + plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + hold off + xlim([1 options_.irf]); + title(deblank(mylist((fig-1)*nstar+plt,:)),'Interpreter','none'); + end + eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig)]); + saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.fig']); + end + if TeX + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for j = 1:nstar + fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist((fig-1)*nstar+j,:)),deblank(mylistTeX((fig-1)*nstar+j,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(fig)); + if options_.relative_irf + fprintf(fidTeX,['\\caption{Relative impulse response' ... + ' functions (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:))); + else + fprintf(fidTeX,['\\caption{Impulse response functions' ... + ' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:))); + end + fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%s}\n',deblank(tit(i,:)),int2str(fig)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end + % close(hh); + end + hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ' figure ' int2str(nbplt) '.']); + m = 0; + for plt = 1:number_of_plots_to_draw-(nbplt-1)*nstar; + m = m+1; + subplot(lr,lc,m); + plot(1:options_.irf,transpose(irfs((nbplt-1)*nstar+plt,:)),'-k','linewidth',1); + hold on + plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + hold off + xlim([1 options_.irf]); + title(deblank(mylist((nbplt-1)*nstar+plt,:)),'Interpreter','none'); + end + eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt)]); + saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt) '.fig']); + end + if TeX + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for j = 1:m + fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist((nbplt-1)*nstar+j,:)),deblank(mylistTeX((nbplt-1)*nstar+j,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(nbplt)); + if options_.relative_irf + fprintf(fidTeX,['\\caption{Relative impulse response functions' ... + ' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:))); + else + fprintf(fidTeX,['\\caption{Impulse response functions' ... + ' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:))); + end + fprintf(fidTeX,'\\label{Fig:IRF:%s:%s}\n',deblank(tit(i,:)),int2str(nbplt)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end + % close(hh); + end + end + end + iter_ = olditer; + if TeX + fprintf(fidTeX,' \n'); + fprintf(fidTeX,'%% End Of TeX file. \n'); + fclose(fidTeX); + end + end + end + + if options_.SpectralDensity == 1 + [omega,f] = UnivariateSpectralDensity(oo_.dr,var_list); + end + + +options_ = options_old; diff --git a/mex/2007b/simulate.mexw32 b/mex/2007b/simulate.mexw32 index d8953dfcd..466daec81 100755 Binary files a/mex/2007b/simulate.mexw32 and b/mex/2007b/simulate.mexw32 differ diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc index dc9aa7545..43cb9e4ec 100644 --- a/mex/sources/simulate/Interpreter.cc +++ b/mex/sources/simulate/Interpreter.cc @@ -193,12 +193,12 @@ Interpreter::compute_block_time() /*throw(EvalException)*/ mexEvalString("drawnow;"); #endif y[(it_+lag)*y_size+var] = Stack.top(); +#ifdef DEBUGC if(var==153) { mexPrintf(" FSTP y[var=%d,time=%d,lag=%d,%d]=%f\n",var,it_,lag,(it_+lag)*y_size+var,y[(it_+lag)*y_size+var]); mexEvalString("drawnow;"); } -#ifdef DEBUGC mexPrintf("%f\n",y[(it_+lag)*y_size+var]); mexEvalString("drawnow;"); #endif diff --git a/mex/sources/simulate/SparseMatrix.cc b/mex/sources/simulate/SparseMatrix.cc index a0f91337d..69b60e2fe 100644 --- a/mex/sources/simulate/SparseMatrix.cc +++ b/mex/sources/simulate/SparseMatrix.cc @@ -601,6 +601,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map< mexEvalString("drawnow;"); #endif tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]]; + //mexPrintf(" u[%d](%f)*y[%d](%f)=%f",it4->second+u_count_init*t,u[it4->second+u_count_init*t],index_vara[var+Size*(y_kmin+t)],y[index_vara[var+Size*(y_kmin+t)]],u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]]); } } else @@ -610,6 +611,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map< #endif b[eq]=it4->second+u_count_init*t; u[b[eq]]+=tmp_b; + //mexPrintf("u[%d]=%f corr=%f\n",b[eq],u[b[eq]],tmp_b); #ifdef PRINT_OUT mexPrintf("=> b[%d]=%f\n", eq, u[b[eq]]); mexEvalString("drawnow;"); @@ -626,6 +628,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map< mexPrintf("end of Init\n"); mexEvalString("drawnow;"); #endif + //mexEvalString("Init"); mxFree(temp_NZE_R); mxFree(temp_NZE_C); } @@ -641,7 +644,7 @@ void SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, std: #ifdef PRINT_OUT mexPrintf("t=%d\n",t); #endif - int ti_y_kmin=-min( t , y_kmin); + int ti_y_kmin=-min( t , y_kmin); int ti_y_kmax= min( periods-(t+1), y_kmax); it4=IM.begin(); eq=-1; @@ -681,6 +684,7 @@ void SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, std: #endif b[eq]=it4->second+u_count_init*t; u[b[eq]]+=tmp_b; + //mexPrintf("u[%d]=%f\n",b[eq],u[b[eq]]); #ifdef PRINT_OUT mexPrintf("=> b[%d]=%f\n", eq, u[b[eq]]); #endif @@ -691,6 +695,7 @@ void SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, std: it4++; } } + //mexPrintf("ShortInit\n"); } diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index b21c48c1f..205122d70 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -31,158 +31,346 @@ using namespace std; #include "BlockTriangular.hh" //------------------------------------------------------------------------------ +/*BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) : + symbol_table(symbol_table_arg), + normalization(symbol_table_arg)*/ BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) : symbol_table(symbol_table_arg), - normalization(symbol_table_arg) + normalization(symbol_table_arg), + incidencematrix(symbol_table_arg) { bt_verbose = 0; ModelBlock = NULL; - Model_Max_Lead = 0; - Model_Max_Lag = 0; periods = 0; } + +IncidenceMatrix::IncidenceMatrix(const SymbolTable &symbol_table_arg) : + symbol_table(symbol_table_arg) +{ + Model_Max_Lead = Model_Max_Lead_Endo = Model_Max_Lead_Exo = 0; + Model_Max_Lag = Model_Max_Lag_Endo = Model_Max_Lag_Exo = 0; + +} //------------------------------------------------------------------------------ //For a lead or a lag build the Incidence Matrix structures List_IM* -BlockTriangular::Build_IM(int lead_lag) +IncidenceMatrix::Build_IM(int lead_lag, SymbolType type) { - List_IM* pIM = new List_IM; int i; - Last_IM->pNext = pIM; - pIM->IM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(pIM->IM[0])); - for(i = 0;i < endo_nbr*endo_nbr;i++) - pIM->IM[i] = 0; - pIM->lead_lag = lead_lag; - if(lead_lag > 0) + List_IM* pIM = new List_IM; + if(type==eEndogenous) { - if(lead_lag > Model_Max_Lead) - Model_Max_Lead = lead_lag; + Last_IM->pNext = pIM; + pIM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) + pIM->IM[i] = 0; + pIM->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Endo) + { + Model_Max_Lead_Endo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Endo) + { + Model_Max_Lag_Endo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM->pNext = NULL; + Last_IM = pIM; } else - { - if( -lead_lag > Model_Max_Lag) - Model_Max_Lag = -lead_lag; + { //eExogenous + Last_IM_X->pNext = pIM; + pIM->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0])); + for(i = 0;i < symbol_table.exo_nbr*symbol_table.endo_nbr;i++) + pIM->IM[i] = 0; + pIM->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Exo) + { + Model_Max_Lead_Exo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Exo) + { + Model_Max_Lag_Exo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM->pNext = NULL; + Last_IM_X = pIM; } - pIM->pNext = NULL; - Last_IM = pIM; return (pIM); } //------------------------------------------------------------------------------ // initialize all the incidence matrix structures void -BlockTriangular::init_incidence_matrix(int nb_endo) +IncidenceMatrix::init_incidence_matrix() { int i; - endo_nbr = nb_endo; First_IM = new List_IM; - First_IM->IM = (bool*)malloc(nb_endo * nb_endo * sizeof(First_IM->IM[0])); - for(i = 0;i < nb_endo*nb_endo;i++) + First_IM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(First_IM->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) First_IM->IM[i] = 0; First_IM->lead_lag = 0; First_IM->pNext = NULL; Last_IM = First_IM; - //cout << "init_incidence_matrix done \n"; + + First_IM_X = new List_IM; + First_IM_X->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(First_IM_X->IM[0])); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.exo_nbr;i++) + First_IM_X->IM[i] = 0; + First_IM_X->lead_lag = 0; + First_IM_X->pNext = NULL; + Last_IM_X = First_IM_X; + } void -BlockTriangular::Free_IM(List_IM* First_IM) const +IncidenceMatrix::Free_IM() const { -#ifdef DEBUG - cout << "Free_IM\n"; -#endif List_IM *Cur_IM, *SFirst_IM; Cur_IM = SFirst_IM = First_IM; while(Cur_IM) { - First_IM = Cur_IM->pNext; + SFirst_IM = Cur_IM->pNext; free(Cur_IM->IM); delete Cur_IM; - Cur_IM = First_IM; + Cur_IM = SFirst_IM; + } + Cur_IM = SFirst_IM = First_IM_X; + while(Cur_IM) + { + SFirst_IM = Cur_IM->pNext; + free(Cur_IM->IM); + delete Cur_IM; + Cur_IM = SFirst_IM; } - //free(SFirst_IM); - //delete SFirst_IM; } //------------------------------------------------------------------------------ -// Return the inceidence matrix related to a lead or a lag +// Return the incidence matrix related to a lead or a lag List_IM* -BlockTriangular::Get_IM(int lead_lag) +IncidenceMatrix::Get_IM(int lead_lag, SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type==eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag)) Cur_IM = Cur_IM->pNext; return (Cur_IM); } bool* -BlockTriangular::bGet_IM(int lead_lag) const +IncidenceMatrix::bGet_IM(int lead_lag, SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type==eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag)) { Cur_IM = Cur_IM->pNext; } - if((Cur_IM->lead_lag != lead_lag) || (Cur_IM==NULL)) - { - cout << "the incidence matrix with lag " << lead_lag << " does not exist !!"; - exit(EXIT_FAILURE); - } - return (Cur_IM->IM); + if(Cur_IM) + return (Cur_IM->IM); + else + return(0); } //------------------------------------------------------------------------------ // Fill the incidence matrix related to a lead or a lag void -BlockTriangular::fill_IM(int equation, int variable_endo, int lead_lag) +IncidenceMatrix::fill_IM(int equation, int variable, int lead_lag, SymbolType type) { List_IM* Cur_IM; - //cout << "equation=" << equation << " variable_endo=" << variable_endo << " lead_lag=" << lead_lag << "\n"; - Cur_IM = Get_IM(lead_lag); - if(equation >= endo_nbr) + Cur_IM = Get_IM(lead_lag, type); + if(equation >= symbol_table.endo_nbr) { - cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; - system("PAUSE"); + cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << symbol_table.endo_nbr << ")\n"; exit(EXIT_FAILURE); } if (!Cur_IM) - Cur_IM = Build_IM(lead_lag); - Cur_IM->IM[equation*endo_nbr + variable_endo] = 1; + Cur_IM = Build_IM(lead_lag, type); + if(type==eEndogenous) + Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 1; + else + Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 1; } //------------------------------------------------------------------------------ // unFill the incidence matrix related to a lead or a lag void -BlockTriangular::unfill_IM(int equation, int variable_endo, int lead_lag) +IncidenceMatrix::unfill_IM(int equation, int variable, int lead_lag, SymbolType type) { List_IM* Cur_IM; - //cout << "lead_lag=" << lead_lag << "\n"; - Cur_IM = Get_IM(lead_lag); - /*if(equation >= endo_nbr) - { - cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; - system("PAUSE"); - exit(EXIT_FAILURE); - }*/ + Cur_IM = Get_IM(lead_lag, type); if (!Cur_IM) - Cur_IM = Build_IM(lead_lag); - Cur_IM->IM[equation*endo_nbr + variable_endo] = 0; - /*system("pause");*/ + Cur_IM = Build_IM(lead_lag, type); + if(type==eEndogenous) + Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 0; + else + Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 0; +} + +List_IM* +IncidenceMatrix::Get_First(SymbolType type) const +{ + if(type==eEndogenous) + return(First_IM); + else + return(First_IM_X); } +//------------------------------------------------------------------------------ +//For a lead or a lag build the Incidence Matrix structures +/*List_IM* +IncidenceMatrix::Build_IM_X(int lead_lag) +{ + List_IM* pIM_X = new List_IM; + int i; + Last_IM_X->pNext = pIM_X; + pIM_X->IM = (bool*)malloc(exo_nbr * endo_nbr * sizeof(pIM_X->IM[0])); + for(i = 0;i < exo_nbr*endo_nbr;i++) + pIM_X->IM[i] = 0; + pIM_X->lead_lag = lead_lag; + if(lead_lag > 0) + { + if(lead_lag > Model_Max_Lead_Exo) + { + Model_Max_Lead_Exo = lead_lag; + if(lead_lag > Model_Max_Lead) + Model_Max_Lead = lead_lag; + } + } + else + { + if( -lead_lag > Model_Max_Lag_Exo) + { + Model_Max_Lag_Exo = -lead_lag; + if(-lead_lag > Model_Max_Lag) + Model_Max_Lag = -lead_lag; + } + } + pIM_X->pNext = NULL; + Last_IM_X = pIM_X; + return (pIM_X); +} + +//------------------------------------------------------------------------------ +// initialize all the incidence matrix structures +void +IncidenceMatrix::init_incidence_matrix_X(int nb_exo) +{ + int i; + //cout << "init_incidence_matrix_X nb_exo = " << nb_exo << " endo_nbr=" << endo_nbr << "\n"; + exo_nbr = nb_exo; + First_IM_X = new List_IM; + First_IM_X->IM = (bool*)malloc(nb_exo * endo_nbr * sizeof(First_IM_X->IM[0])); + for(i = 0;i < endo_nbr*nb_exo;i++) + First_IM_X->IM[i] = 0; + First_IM_X->lead_lag = 0; + First_IM_X->pNext = NULL; + Last_IM_X = First_IM_X; +} + + +void +IncidenceMatrix::Free_IM_X(List_IM* First_IM_X) const +{ + List_IM *Cur_IM_X, *SFirst_IM_X; + Cur_IM_X = SFirst_IM_X = First_IM_X; + while(Cur_IM_X) + { + First_IM_X = Cur_IM_X->pNext; + free(Cur_IM_X->IM); + delete Cur_IM_X; + Cur_IM_X = First_IM_X; + } +} + + +//------------------------------------------------------------------------------ +// Return the inceidence matrix related to a lead or a lag +List_IM* +IncidenceMatrix::Get_IM_X(int lead_lag) +{ + List_IM* Cur_IM_X; + Cur_IM_X = First_IM_X; + while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag)) + Cur_IM_X = Cur_IM_X->pNext; + return (Cur_IM_X); +} + +bool* +IncidenceMatrix::bGet_IM_X(int lead_lag) const +{ + List_IM* Cur_IM_X; + Cur_IM_X = First_IM_X; + while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag)) + { + Cur_IM_X = Cur_IM_X->pNext; + } + if(Cur_IM_X) + return (Cur_IM_X->IM); + else + return(0); +} + + +//------------------------------------------------------------------------------ +// Fill the incidence matrix related to a lead or a lag +void +IncidenceMatrix::fill_IM_X(int equation, int variable_exo, int lead_lag) +{ + List_IM* Cur_IM_X; + Cur_IM_X = Get_IM_X(lead_lag); + if(equation >= endo_nbr) + { + cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n"; + exit(EXIT_FAILURE); + } + if (!Cur_IM_X) + { + Cur_IM_X = Build_IM_X(lead_lag); + } + Cur_IM_X->IM[equation*exo_nbr + variable_exo] = true; +} +*/ + //------------------------------------------------------------------------------ //Print azn incidence matrix void -BlockTriangular::Print_SIM(bool* IM, int n) const +IncidenceMatrix::Print_SIM(bool* IM, SymbolType type) const { - int i, j; - for(i = 0;i < n;i++) + int i, j, n; + if(type == eEndogenous) + n = symbol_table.endo_nbr; + else + n = symbol_table.exo_nbr; + for(i = 0;i < symbol_table.endo_nbr;i++) { cout << " "; for(j = 0;j < n;j++) @@ -194,15 +382,18 @@ BlockTriangular::Print_SIM(bool* IM, int n) const //------------------------------------------------------------------------------ //Print all incidence matrix void -BlockTriangular::Print_IM(int n) const +IncidenceMatrix::Print_IM(SymbolType type) const { List_IM* Cur_IM; - Cur_IM = First_IM; + if(type == eEndogenous) + Cur_IM = First_IM; + else + Cur_IM = First_IM_X; cout << "-------------------------------------------------------------------\n"; while(Cur_IM) { cout << "Incidence matrix for lead_lag = " << Cur_IM->lead_lag << "\n"; - Print_SIM(Cur_IM->IM, n); + Print_SIM(Cur_IM->IM, type); Cur_IM = Cur_IM->pNext; } } @@ -211,7 +402,7 @@ BlockTriangular::Print_IM(int n) const //------------------------------------------------------------------------------ // Swap rows and columns of the incidence matrix void -BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) +IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const { int tmp_i, j; bool tmp_b; @@ -269,15 +460,13 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n if ((k == 1) && IM0[Index_Equ_IM[i].index*n + Index_Var_IM[l].index]) { modifie = 1; - swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n); + incidencematrix.swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n); (*prologue)++; } } } *epilogue = 0; modifie = 1; - /*print_SIM(IM,n); - print_SIM(IM*/ while(modifie) { modifie = 0; @@ -295,7 +484,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n if ((k == 1) && IM0[Index_Equ_IM[l].index*n + Index_Var_IM[i].index]) { modifie = 1; - swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n); + incidencematrix.swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n); (*epilogue)++; } } @@ -306,11 +495,12 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n void BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock) { - int i, j, k, l, ls, m, i_1, Lead, Lag, size_list_lead_var, first_count_equ, i1; - int *list_lead_var, *tmp_size, *tmp_var, *tmp_endo, nb_lead_lag_endo; + int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1; + int *tmp_size, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_exo, tmp_nb_exo, nb_lead_lag_endo; List_IM *Cur_IM; bool *IM, OK; ModelBlock->Periods = periods; + int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo; if ((type == PROLOGUE) || (type == EPILOGUE)) { for(i = 0;i < size;i++) @@ -321,128 +511,221 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN; ModelBlock->Block_List[*count_Block].Temporary_terms=new temporary_terms_type (); ModelBlock->Block_List[*count_Block].Temporary_terms->clear(); - list_lead_var = (int*)malloc(Model_Max_Lead * endo_nbr * sizeof(int)); - size_list_lead_var = 0; - tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); + tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); tmp_var = (int*)malloc(sizeof(int)); - memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); - memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); + tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); nb_lead_lag_endo = Lead = Lag = 0; - Cur_IM = First_IM; + Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + Cur_IM = incidencematrix.Get_First(eEndogenous); while(Cur_IM) { k = Cur_IM->lead_lag; - i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; if(k > 0) { - if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) { nb_lead_lag_endo++; - tmp_size[Model_Max_Lag + k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if(k > Lead) - { - Lead = k; - list_lead_var[size_list_lead_var] = Index_Var_IM[*count_Equ].index + size * (k - 1); - size_list_lead_var++; - } + Lead = k; } } else { k = -k; - if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index]) + if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index]) { - tmp_size[Model_Max_Lag - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; nb_lead_lag_endo++; if(k > Lag) - { - Lag = k; - } + Lag = k; } } Cur_IM = Cur_IM->pNext; } + + Lag_Endo = Lag; + Lead_Endo = Lead; + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + + + Cur_IM = incidencematrix.Get_First(eExogenous); + k = Cur_IM->lead_lag; + while(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j]) + { + if(!tmp_exo[j]) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + if(k>0 && k>Lead_Exo) + Lead_Exo = k; + else if(k<0 && (-k)>Lag_Exo) + Lag_Exo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++; + } + Cur_IM = Cur_IM->pNext; + if(Cur_IM) + k = Cur_IM->lead_lag; + } + + + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } + + ModelBlock->Block_List[*count_Block].nb_exo_det = 0; + ModelBlock->Block_List[*count_Block].Max_Lag = Lag; ModelBlock->Block_List[*count_Block].Max_Lead = Lead; - free(list_lead_var); + ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo; + ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo; + ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo; + ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo; ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(sizeof(int)); - ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Equation[0] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[0] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].Variable_Sorted[0] = -1; - ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = 0; if ((Lead > 0) && (Lag > 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE; else if((Lead > 0) && (Lag == 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_BACKWARD_SIMPLE; else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; + + Cur_IM = incidencematrix.Get_First(eExogenous); + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + while(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j] && (!tmp_exo[j])) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + Cur_IM = Cur_IM->pNext; + } + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; - if(nb_lead_lag_endo) - { - ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - } + + + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } ls = l = 1; i1 = 0; for(int li = 0;li < Lead + Lag + 1;li++) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[Model_Max_Lag - Lag + li]; - if(tmp_size[Model_Max_Lag - Lag + li]) + if(incidencematrix.Model_Max_Lag_Endo - Lag + li>=0) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[Model_Max_Lag - Lag + li]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int)); + //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].size=" << tmp_size[Model_Max_Lag_Endo - Lag + li] << "\n"; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l; - IM = bGet_IM(li - Lag); - if(IM == NULL) + IM = incidencematrix.bGet_IM(li - Lag, eEndogenous); + if(IM) { - cout << "Error IM(" << li - Lag << ") doesn't exist\n"; - exit(EXIT_FAILURE); - } - if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr] && nb_lead_lag_endo) - { - ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = li - Lag; - tmp_var[0] = i1; - i1++; - } - m = 0; - i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr; - if(IM[Index_Var_IM[*count_Equ].index + i_1]) - { - if(li == Lag) + if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*symbol_table.endo_nbr] && nb_lead_lag_endo) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls; - ls++; + tmp_var[0] = i1; + i1++; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]]; - l++; - m++; + m = 0; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; + if(IM[Index_Var_IM[*count_Equ].index + i_1]) + { + if(li == Lag) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls; + ls++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; + l++; + m++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; - } + } + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = 0; + if(incidencematrix.Model_Max_Lag_Exo - Lag + li>=0) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li]; + //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].Exogenous= malloc(" << tmp_size_exo[Model_Max_Lag_Exo - Lag + li] << ")\n"; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int)); + IM = incidencematrix.bGet_IM(li - Lag, eExogenous); + if(IM) + { + m = 0; + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr; + for(k = 0; kBlock_List[*count_Block].Exogenous[k]+i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous[m] = k; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X[m] = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index[m] = Index_Equ_IM[*count_Equ].index; + m++; + } + } + } + } + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = 0; } (*count_Equ)++; (*count_Block)++; free(tmp_size); + free(tmp_size_exo); free(tmp_endo); + free(tmp_exo); free(tmp_var); } } @@ -456,25 +739,25 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN; ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); - ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); + //ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int)); Lead = Lag = 0; first_count_equ = *count_Equ; tmp_var = (int*)malloc(size * sizeof(int)); - tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int)); - memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); - memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int)); + tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); nb_lead_lag_endo = 0; + Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + for(i = 0;i < size;i++) { ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[i] = Index_Var_IM[*count_Equ].index; - ModelBlock->Block_List[*count_Block].Variable_Sorted[i] = -1; - ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block; - ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = i; - Cur_IM = First_IM; + Cur_IM = incidencematrix.Get_First(eEndogenous); i_1 = Index_Var_IM[*count_Equ].index; while(Cur_IM) { @@ -484,12 +767,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { for(j = 0;j < size;j++) { - if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr]) + if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { - tmp_size[Model_Max_Lag + k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if (!OK) { - tmp_endo[Model_Max_Lag + k]++; + tmp_endo[incidencematrix.Model_Max_Lag + k]++; nb_lead_lag_endo++; OK = true; } @@ -503,12 +786,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc k = -k; for(j = 0;j < size;j++) { - if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr]) + if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { - tmp_size[Model_Max_Lag - k]++; + tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++; if (!OK) { - tmp_endo[Model_Max_Lag - k]++; + tmp_endo[incidencematrix.Model_Max_Lag - k]++; nb_lead_lag_endo++; OK = true; } @@ -537,87 +820,158 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; } + + + Lag_Endo = Lag; + Lead_Endo = Lead; + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + tmp_nb_exo = 0; + for(i = 0;i < size;i++) + { + Cur_IM = incidencematrix.Get_First(eExogenous); + k = Cur_IM->lead_lag; + while(Cur_IM) + { + i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.exo_nbr; + for(j=0;jIM[i_1 + j]) + { + if(!tmp_exo[j]) + { + tmp_exo[j] = 1; + tmp_nb_exo++; + } + if(k>0 && k>Lead_Exo) + Lead_Exo = k; + else if(k<0 && (-k)>Lag_Exo) + Lag_Exo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++; + } + Cur_IM = Cur_IM->pNext; + if(Cur_IM) + k = Cur_IM->lead_lag; + } + } + + + ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo; + ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int)); + k = 0; + for(j=0;jBlock_List[*count_Block].Exogenous[k] = j; + k++; + } + + ModelBlock->Block_List[*count_Block].nb_exo_det = 0; + ModelBlock->Block_List[*count_Block].Max_Lag = Lag; ModelBlock->Block_List[*count_Block].Max_Lead = Lead; + ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo; + ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo; + ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo; + ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo; ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); ls = l = size; i1 = 0; ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; - ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int)); - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int)); for(i = 0;i < Lead + Lag + 1;i++) { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[Model_Max_Lag - Lag + i]; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int)); - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; - IM = bGet_IM(i - Lag); - if(IM != NULL) + if(incidencematrix.Model_Max_Lag_Endo-Lag+i>=0) { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + 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)); } else + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = 0; + if(incidencematrix.Model_Max_Lag_Exo-Lag+i>=0) { - cout << "Error IM(" << i - Lag << ") doesn't exist\n"; - exit(EXIT_FAILURE); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int)); } - for(j = first_count_equ;j < size + first_count_equ;j++) + else + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = 0; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; + IM = incidencematrix.bGet_IM(i - Lag, eEndogenous); + if(IM) { - i_1 = Index_Var_IM[j].index; - m = 0; - for(k = first_count_equ;k < size + first_count_equ;k++) - if(IM[i_1 + Index_Equ_IM[k].index*endo_nbr]) - m++; - if(m > 0) + for(j = first_count_equ;j < size + first_count_equ;j++) { - ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = i_1; - ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = i - Lag; - tmp_var[j - first_count_equ] = i1; - i1++; + i_1 = Index_Var_IM[j].index; + m = 0; + for(k = first_count_equ;k < size + first_count_equ;k++) + if(IM[i_1 + Index_Equ_IM[k].index*symbol_table.endo_nbr]) + m++; + if(m > 0) + { + tmp_var[j - first_count_equ] = i1; + i1++; + } + } + m = 0; + for(j = first_count_equ;j < size + first_count_equ;j++) + { + i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr; + for(k = first_count_equ;k < size + first_count_equ;k++) + if(IM[Index_Var_IM[k].index + i_1]) + { + if(i == Lag) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; + ls++; + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index; + l++; + m++; + } + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; + } + IM = incidencematrix.bGet_IM(i - Lag, eExogenous); + if(IM) + { + m = 0; + for(j = first_count_equ;j < size + first_count_equ;j++) + { + i_1 = Index_Equ_IM[j].index * symbol_table.exo_nbr; + for(k = 0; kBlock_List[*count_Block].Exogenous[k]+i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous[m] = k; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X[m] = j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j].index; + m++; + } + } } } - m = 0; - for(j = first_count_equ;j < size + first_count_equ;j++) - { - i_1 = Index_Equ_IM[j].index * endo_nbr; - for(k = first_count_equ;k < size + first_count_equ;k++) - if(IM[Index_Var_IM[k].index + i_1]) - { - if(i == Lag) - { - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; - ls++; -#ifdef DEBUG - printf("j=%d ModelBlock->Block_List[%d].Variable[%d]=%d Index_Var_IM[%d].index=%d", j, *count_Block, j - first_count_equ, ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ], k, Index_Var_IM[k].index); - if(ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ] == Index_Var_IM[k].index) - { - ModelBlock->Block_List[*count_Block].Own_Derivative[j - first_count_equ]=l; - printf(" l=%d OK\n",l); - } - else - printf("\n"); -#endif - } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index; - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[k - first_count_equ]]; - l++; - m++; - } - } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; } (*count_Block)++; free(tmp_size); + free(tmp_size_exo); free(tmp_endo); + free(tmp_exo); free(tmp_var); } } @@ -627,11 +981,6 @@ void BlockTriangular::Free_Block(Model_Block* ModelBlock) const { int blk, i; -#ifdef DEBUG - - cout << "Free_Block\n"; -#endif - for(blk = 0;blk < ModelBlock->Size;blk++) { @@ -639,22 +988,26 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const { free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); - free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Exogenous); free(ModelBlock->Block_List[blk].Own_Derivative); - if(ModelBlock->Block_List[blk].Nb_Lead_Lag_Endo) - { - free(ModelBlock->Block_List[blk].variable_dyn_index); - free(ModelBlock->Block_List[blk].variable_dyn_leadlag); - } for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { - free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index); + if(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); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index); + } + if(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); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X); + } } free(ModelBlock->Block_List[blk].IM_lead_lag); delete(ModelBlock->Block_List[blk].Temporary_terms); @@ -663,28 +1016,31 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const { free(ModelBlock->Block_List[blk].Equation); free(ModelBlock->Block_List[blk].Variable); - free(ModelBlock->Block_List[blk].Variable_Sorted); + free(ModelBlock->Block_List[blk].Exogenous); free(ModelBlock->Block_List[blk].Own_Derivative); - free(ModelBlock->Block_List[blk].variable_dyn_index); - free(ModelBlock->Block_List[blk].variable_dyn_leadlag); for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++) { - free(ModelBlock->Block_List[blk].IM_lead_lag[i].u); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].us); - free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); - 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].Var_dyn_Index); + 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); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ); + 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); + } + 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); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index); + free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X); + } } free(ModelBlock->Block_List[blk].IM_lead_lag); delete(ModelBlock->Block_List[blk].Temporary_terms); } } - free(ModelBlock->in_Block_Equ); - free(ModelBlock->in_Block_Var); - free(ModelBlock->in_Equ_of_Block); - free(ModelBlock->in_Var_of_Block); free(ModelBlock->Block_List); free(ModelBlock); free(Index_Equ_IM); @@ -700,40 +1056,19 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int i, j, Nb_TotalBlocks, Nb_RecursBlocks; int count_Block, count_Equ; block_result_t* res; - //List_IM * p_First_IM, *p_Cur_IM, *Cur_IM; Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set)); bool* SIM0, *SIM00; - /*p_First_IM = (List_IM*)malloc(sizeof(*p_First_IM)); - p_Cur_IM = p_First_IM; - Cur_IM = First_IM; - i = endo_nbr * endo_nbr; - while(Cur_IM) - { - p_Cur_IM->lead_lag = Cur_IM->lead_lag; - p_Cur_IM->IM = (bool*)malloc(i * sizeof(bool)); - memcpy ( p_Cur_IM->IM, Cur_IM->IM, i); - Cur_IM = Cur_IM->pNext; - if(Cur_IM) - { - p_Cur_IM->pNext = (List_IM*)malloc(sizeof(*p_Cur_IM)); - p_Cur_IM = p_Cur_IM->pNext; - } - else - p_Cur_IM->pNext = NULL; - }*/ SIM0 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; memcpy(SIM0,IM_0,n*n*sizeof(bool)); Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0); - //cout << "free SIM0=" << SIM0 << "\n"; free(SIM0); if(bt_verbose) { cout << "prologue : " << *prologue << " epilogue : " << *epilogue << "\n"; cout << "IM_0\n"; - Print_SIM(IM_0, n); + incidencematrix.Print_SIM(IM_0, eEndogenous); cout << "IM\n"; - Print_SIM(IM, n); + incidencematrix.Print_SIM(IM, eEndogenous); for(i = 0;i < n;i++) cout << "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index << " Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index << "\n"; } @@ -745,11 +1080,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, if(mixing) { double* max_val=(double*)malloc(n*sizeof(double)); - //cout << "n=" << n << "\n"; memset(max_val,0,n*sizeof(double)); for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) { - //cout << "iter->first.first=" << iter->first.first << "\n"; if(fabs(iter->second)>max_val[iter->first.first]) max_val[iter->first.first]=fabs(iter->second); } @@ -763,17 +1096,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, { int suppress=0; SIM0 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n"; memset(SIM0,0,n*n*sizeof(bool)); SIM00 = (bool*)malloc(n * n * sizeof(bool)); - //cout << "Allocate SIM00=" << SIM00 << " size=" << n * n * sizeof(bool) << "\n"; memset(SIM00,0,n*n*sizeof(bool)); - //cout << "n*n=" << n*n << "\n"; for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ ) { if(fabs(iter->second)>bi) { - //cout << "iter->first.first*n+iter->first.second=" << iter->first.first*n+iter->first.second << "\n"; SIM0[iter->first.first*n+iter->first.second]=1; if(!IM_0[iter->first.first*n+iter->first.second]) { @@ -783,14 +1112,11 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, else suppress++; } - //cout << "n*n=" << n*n << "\n"; for(i = 0;i < n;i++) for(j = 0;j < n;j++) { - //cout << "Index_Equ_IM[i].index * n + Index_Var_IM[j].index=" << Index_Equ_IM[i].index * n + Index_Var_IM[j].index << "\n"; SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index]; } - //cout << "free SIM0=" << SIM0 << "\n"; free(SIM0); if(suppress!=suppressed) { @@ -850,15 +1176,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, cout << " the largest simultaneous block has " << j << " equation(s). \n"; ModelBlock->Size = Nb_TotalBlocks; ModelBlock->Periods = periods; - ModelBlock->in_Block_Equ = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Block_Var = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Equ_of_Block = (int*)malloc(n * sizeof(int)); - ModelBlock->in_Var_of_Block = (int*)malloc(n * sizeof(int)); ModelBlock->Block_List = (Block*)malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks); - blocks.block_result_to_IM(res, IM, *prologue, endo_nbr, Index_Equ_IM, Index_Var_IM); - //Free_IM(p_First_IM); + blocks.block_result_to_IM(res, IM, *prologue, symbol_table.endo_nbr, Index_Equ_IM, Index_Var_IM); count_Equ = count_Block = 0; - //Print_IM(endo_nbr); if (*prologue) Allocate_Block(*prologue, &count_Equ, &count_Block, PROLOGUE, ModelBlock); for(j = 0;j < res->n_sets;j++) @@ -887,11 +1207,11 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ List_IM* Cur_IM; int i; //First create a static model incidence matrix - SIM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM)); - for(i = 0;i < endo_nbr*endo_nbr;i++) + SIM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM)); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) { SIM[i] = 0; - Cur_IM = First_IM; + Cur_IM = incidencematrix.Get_First(eEndogenous); while(Cur_IM) { SIM[i] = (SIM[i]) || (Cur_IM->IM[i]); @@ -901,26 +1221,26 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_ if(bt_verbose) { cout << "incidence matrix for the static model (unsorted) \n"; - Print_SIM(SIM, endo_nbr); + incidencematrix.Print_SIM(SIM, eEndogenous); } - Index_Equ_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Equ_IM)); - for(i = 0;i < endo_nbr;i++) + Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Equ_IM)); + for(i = 0;i < symbol_table.endo_nbr;i++) { Index_Equ_IM[i].index = i; } - Index_Var_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Var_IM)); - for(i = 0;i < endo_nbr;i++) + Index_Var_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Var_IM)); + for(i = 0;i < symbol_table.endo_nbr;i++) { Index_Var_IM[i].index = i; } if(ModelBlock != NULL) Free_Block(ModelBlock); ModelBlock = (Model_Block*)malloc(sizeof(*ModelBlock)); - Cur_IM = Get_IM(0); - SIM_0 = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM_0)); - for(i = 0;i < endo_nbr*endo_nbr;i++) + Cur_IM = incidencematrix.Get_IM(0, eEndogenous); + SIM_0 = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM_0)); + for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++) SIM_0[i] = Cur_IM->IM[i]; - Normalize_and_BlockDecompose(SIM, ModelBlock, endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m); + Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m); free(SIM_0); free(SIM); } diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index c886f83b0..a56f953e7 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -169,7 +169,11 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const { options_list.writeOutput(output); symbol_list.writeOutput("var_list_", output); - output << "info = stoch_simul(var_list_);\n"; + output << "if(options_.model_mode)\n"; + output << " info = stoch_simul_sparse(var_list_);\n"; + output << "else\n"; + output << " info = stoch_simul(var_list_);\n"; + output << "end\n"; } ForecastStatement::ForecastStatement(const SymbolList &symbol_list_arg, diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 592b5d213..8b4bf6163 100644 --- a/preprocessor/ExprNode.cc +++ b/preprocessor/ExprNode.cc @@ -164,6 +164,12 @@ NumConstNode::collectEndogenous(set > &result) const { } +void +NumConstNode::collectExogenous(set > &result) const +{ +} + + VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg) : ExprNode(datatree_arg), symb_id(symb_id_arg), @@ -473,6 +479,14 @@ VariableNode::collectEndogenous(set > &result) const result.insert(make_pair(symb_id, lag)); } +void +VariableNode::collectExogenous(set > &result) const +{ + if (type == eExogenous) + result.insert(make_pair(symb_id, lag)); +} + + UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) : ExprNode(datatree_arg), arg(arg_arg), @@ -864,6 +878,13 @@ UnaryOpNode::collectEndogenous(set > &result) const arg->collectEndogenous(result); } +void +UnaryOpNode::collectExogenous(set > &result) const +{ + arg->collectExogenous(result); +} + + BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, BinaryOpcode op_code_arg, const NodeID arg2_arg) : ExprNode(datatree_arg), @@ -1322,6 +1343,14 @@ BinaryOpNode::collectEndogenous(set > &result) const arg2->collectEndogenous(result); } + +void +BinaryOpNode::collectExogenous(set > &result) const +{ + arg1->collectExogenous(result); + arg2->collectExogenous(result); +} + TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) : ExprNode(datatree_arg), @@ -1334,7 +1363,7 @@ TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg, // Non-null derivatives are the union of those of the arguments // Compute set union of arg{1,2,3}->non_null_derivatives - set non_null_derivatives_tmp; + set non_null_derivatives_tmp; set_union(arg1->non_null_derivatives.begin(), arg1->non_null_derivatives.end(), arg2->non_null_derivatives.begin(), @@ -1590,6 +1619,15 @@ TrinaryOpNode::collectEndogenous(set > &result) const arg3->collectEndogenous(result); } +void +TrinaryOpNode::collectExogenous(set > &result) const +{ + arg1->collectExogenous(result); + arg2->collectExogenous(result); + arg3->collectExogenous(result); +} + + UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg, const vector &arguments_arg) : @@ -1650,6 +1688,15 @@ UnknownFunctionNode::collectEndogenous(set > &result) const (*it)->collectEndogenous(result); } +void +UnknownFunctionNode::collectExogenous(set > &result) const +{ + for(vector::const_iterator it = arguments.begin(); + it != arguments.end(); it++) + (*it)->collectExogenous(result); +} + + double UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (EvalException) { diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index 495b030a2..f0d064250 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -69,7 +69,7 @@ ModFile::checkPass() exit(EXIT_FAILURE); } - if (mod_file_struct.simul_present && stochastic_statement_present) + if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.mode==0) { cerr << "ERROR: A .mod file cannot contain both a simul command and one of {stoch_simul, estimation, forecast, osr, ramsey_policy}" << endl; exit(EXIT_FAILURE); diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 59fda0327..116c61f04 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -53,9 +53,10 @@ ModelTree::equation_number() const void ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, - const temporary_terms_type &temporary_terms) const + const temporary_terms_type &temporary_terms, + SymbolType type) const { - first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag))); + first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(type, symb_id, lag))); if (it != first_derivatives.end()) (it->second)->writeOutput(output, output_type, temporary_terms); else @@ -67,14 +68,9 @@ ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, { first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag))); if (it != first_derivatives.end()) - { - /*NodeID Id = it->second;*/ - (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx); - } + (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx); else - { - code_file.write(&FLDZ, sizeof(FLDZ)); - } + code_file.write(&FLDZ, sizeof(FLDZ)); } @@ -219,6 +215,35 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t } } + +void +ModelTree::BuildIncidenceMatrix() +{ + set > endogenous, exogenous; + for(int eq = 0; eq < (int) equations.size(); eq++) + { + BinaryOpNode *eq_node = equations[eq]; + endogenous.clear(); + NodeID Id = eq_node->arg1; + Id->collectEndogenous(endogenous); + Id = eq_node->arg2; + Id->collectEndogenous(endogenous); + for(set >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, it_endogenous->second, eEndogenous); + } + exogenous.clear(); + Id = eq_node->arg1; + Id->collectExogenous(exogenous); + Id = eq_node->arg2; + Id->collectExogenous(exogenous); + for(set >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + { + block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, it_exogenous->second, eExogenous); + } + } +} + void ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const { @@ -244,7 +269,7 @@ void ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) { map reference_count, first_occurence; - int i, j, m, eq, var, lag/*, prev_size=0*/; + int i, j, m, eq, var, lag; temporary_terms_type vect; ostringstream tmp_output; BinaryOpNode *eq_node; @@ -265,6 +290,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) tmp_output.str(""); lhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms); tmp_s << "y[Per_y_+" << ModelBlock->Block_List[j].Variable[0] << "]"; + //Determine whether the equation could be evaluated rather than to be solved if (tmp_output.str()==tmp_s.str()) { if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE) @@ -285,6 +311,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) } } } + // Compute the temporary terms reordered for(i = 0;i < ModelBlock->Block_List[j].Size;i++) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -335,12 +362,11 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) for(second_derivatives_type::iterator it = second_derivatives.begin(); it != second_derivatives.end(); it++) it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, 0, ModelBlock, map_idx); - /*New*/ + // Add a mapping form node ID to temporary terms order j=0; for(temporary_terms_type::const_iterator it = temporary_terms.begin(); it != temporary_terms.end(); it++) map_idx[(*it)->idx]=j++; - /*EndNew*/ } void @@ -355,6 +381,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock ostringstream Uf[symbol_table.endo_nbr]; map reference_count; int prev_Simulation_Type=-1, count_derivates=0; + int jacobian_max_endo_col; temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); //---------------------------------------------------------------------- //Temporary variables declaration @@ -416,9 +443,9 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n"; else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, g1, g2, g3, y_index, jacobian_eval)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n"; else - output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods, g1, g2, g3)\n"; + output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, periods, jacobian_eval, g1, g2, g3, y_kmin, y_size)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -434,13 +461,13 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock temporary_terms_type tt2; - if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE) + if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE) { int nze; for(nze=0,m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size; - //output << " Jacobian_Size=" << ModelBlock->Block_List[j].Size << "*(y_kmin+" << ModelBlock->Block_List[j].Max_Lead << " +periods);\n"; - //output << " g1=spalloc( y_size*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; + output << " g2=0;g3=0;\n"; + output << " b = [];\n"; output << " for it_ = y_kmin+1:(periods+y_kmin)\n"; output << " Per_y_=it_*y_size;\n"; output << " Per_J_=(it_-y_kmin-1)*y_size;\n"; @@ -467,7 +494,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock // 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; @@ -502,7 +528,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock goto end; case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: - Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")"; + Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")"; output << sps << "residual(" << i+1 << ") = ("; goto end; case SOLVE_TWO_BOUNDARIES_COMPLETE: @@ -543,23 +569,42 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { if(ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]==ModelBlock->Block_List[j].Variable[0]) { - //output << " g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - //output << " g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - output << " g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", " << ModelBlock->Block_List[j].Variable[0]+1 + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms); + output << " g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", " + << ModelBlock->Block_List[j].Variable[0]+1 + + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr + << ")="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) - << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) << ") " << ModelBlock->Block_List[j].Variable[0]+1 << ", equation=" << ModelBlock->Block_List[j].Equation[0]+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; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eq+1 << ", " + << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) { output << " else\n"; output << " g1="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms); + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) << ") " << ModelBlock->Block_List[j].Variable[0]+1 @@ -576,6 +621,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: count_derivates++; + output << " b = [];\n"; for(m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++) { k=m-ModelBlock->Block_List[j].Max_Lag; @@ -583,17 +629,35 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " g1(" << eqr+1 << ", " << varr+1+m*ModelBlock->Block_List[j].Size << ")="; - writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " g1(" << eq+1 << ", " << var+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << /*variable_table.getLag(variable_table.getSymbolID(var))*/k << ") " << var+1 << ", equation=" << eq+1 << endl; } } - output << " else\n"; + 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; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eq+1 << ", " << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + output << " else" << endl; + m=ModelBlock->Block_List[j].Max_Lag; for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { @@ -601,22 +665,20 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock 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]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << ", " << varr+1 << ")*y(it_, " << var+1 << ")"; - //output << " u(" << u+1 << ") = "; output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms); + writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 << ", equation=" << eq+1 << endl; } - output << " end;\n"; for(i = 0;i < ModelBlock->Block_List[j].Size;i++) output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + output << " end;\n"; break; case SOLVE_TWO_BOUNDARIES_SIMPLE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " g2=0;g3=0;\n"; + output << " if ~jacobian_eval" << endl; 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; @@ -624,7 +686,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - //int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[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 (k==0) @@ -637,14 +698,14 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; //output << " u(" << u+1 << "+Per_u_) = "; if(k==0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; else if(k==1) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = "; else if(k>0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = "; else if(k<0) - output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; - writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms); + output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; @@ -656,7 +717,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock } for(i = 0;i < ModelBlock->Block_List[j].Size;i++) { - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; #ifdef CONDITION output << " if (fabs(condition(" << i+1 << "))Block_List[j].Size;i++) output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; #endif + + output << " else" << endl; + 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; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; + 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]; + output << " g1(" << eqr+1 << ", " << varr+1+(m-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size; + 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; + for(i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i]; + int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; + output << " g1(" << eqr+1 << ", " + << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); + output << "; % variable (exogenous)=" << symbol_table.getNameByID(eExogenous, var) + << "(" << k << ") " << var+1 << " " << varr+1 + << ", equation=" << eq+1 << endl; + } + } + output << " end;\n"; output << " end;\n"; break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } @@ -689,7 +789,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock void ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &static_basename) const { - int i,j,k,m, var, eq; + int i,j,k,m, var, eq, g1_index = 1; string tmp_s, sps; ostringstream tmp_output, global_output; NodeID lhs=NULL, rhs=NULL; @@ -734,9 +834,15 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R )) - skip_the_head=true; + { + skip_the_head=true; + g1_index++; + } else - skip_the_head=false; + { + skip_the_head=false; + g1_index = 1; + } if (!skip_the_head) { if (j>0) @@ -749,9 +855,9 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R ) - output << "function [y] = " << static_basename << "_" << j+1 << "(y, x)\n"; + output << "function [y, g1] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n"; else - output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n"; + output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n"; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << " //" << endl @@ -773,14 +879,17 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode memset(IM, 0, n*n*sizeof(bool)); for(m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++) { - IMl=block_triangular.bGet_IM(m); - for(i=0;iBlock_List[j].Equation[i]; - for(k=0;kBlock_List[j].Variable[k]; - IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; + eq=ModelBlock->Block_List[j].Equation[i]; + for(k=0;kBlock_List[j].Variable[k]; + IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var]; + } } } } @@ -814,7 +923,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode // 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); + //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; @@ -862,108 +971,114 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode } } // The Jacobian if we have to solve the block - if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R - && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R) + output << " " << sps << "% Jacobian " << endl; + switch(ModelBlock->Block_List[j].Simulation_Type) { - output << " " << sps << "% Jacobian " << endl; - switch(ModelBlock->Block_List[j].Simulation_Type) + case EVALUATE_BACKWARD: + case EVALUATE_FORWARD: + case EVALUATE_BACKWARD_R: + case EVALUATE_FORWARD_R: + output << " if(jacobian_eval)\n"; + output << " g1( " << g1_index << ", " << g1_index << ")="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) + << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + output << " end\n"; + break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: + output << " g1(1)="; + writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) + << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) + << ") " << ModelBlock->Block_List[j].Variable[0]+1 + << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; + break; + case SOLVE_BACKWARD_COMPLETE: + case SOLVE_FORWARD_COMPLETE: + output << " g2=0;g3=0;\n"; + m=ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: - output << " g1(1)="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) - << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0]+1 - << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; - break; - case SOLVE_BACKWARD_COMPLETE: - case SOLVE_FORWARD_COMPLETE: - output << " g2=0;g3=0;\n"; - m=ModelBlock->Block_List[j].Max_Lag; + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; + 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(" << var+1 << ")"; + output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; + writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; + break; + case SOLVE_TWO_BOUNDARIES_COMPLETE: + case SOLVE_TWO_BOUNDARIES_SIMPLE: + output << " g2=0;g3=0;\n"; + 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; for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) { int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - //int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[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]] << "-u(" << u << ")*y(Per_y_+" << var << ")"; - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")"; - //output << " u(" << u+1 << ") = "; - output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms); + 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]=1; + } + 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) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 + << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; - break; - case SOLVE_TWO_BOUNDARIES_COMPLETE: - output << " g2=0;g3=0;\n"; - 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; - for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - //int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[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]; - //output << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n"; - 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]=1; - } - output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; - writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << k << ") " << var+1 - << ", equation=" << eq+1 << endl; #ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++) - { - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; -#ifdef CONDITION - output << " if (fabs(condition(" << i+1 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) - { - k=m-ModelBlock->Block_List[j].Max_Lag; - for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) - { - int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; - int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; - int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; - int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; - output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; - } - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; -#endif - break; } + output << " if(~jacobian_eval)\n"; + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + { + output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; +#ifdef CONDITION + output << " if (fabs(condition(" << i+1 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) + { + k=m-ModelBlock->Block_List[j].Max_Lag; + for(i=0;iBlock_List[j].IM_lead_lag[m].size;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; + int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i]; + int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i]; + output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n"; + } + } + for(i = 0;i < ModelBlock->Block_List[j].Size;i++) + output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n"; +#endif + break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; free(IM); } output << "return;\n\n\n"; - //free(IM); } @@ -1031,7 +1146,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; } ModelBlock_Aggregated_Count++; - //cout << "ModelBlock_Aggregated_Count=" << ModelBlock_Aggregated_Count << "\n"; //For each block j=0; for(k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++) @@ -1044,7 +1158,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl code_file.write(reinterpret_cast(&v),sizeof(v)); v=ModelBlock->Block_List[j].Simulation_Type; code_file.write(reinterpret_cast(&v),sizeof(v)); - //cout << "FBEGINBLOCK j=" << j << " size=" << ModelBlock_Aggregated_Number[k0] << " type=" << v << "\n"; for(k=0; kBlock_List[j].Size;i++) @@ -1079,7 +1192,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } for(k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++) { - //For a block composed of a single equation determines wether we have to evaluate or to solve the equation + //For a block composed of a single equation determines whether we have to evaluate or to solve the equation if (ModelBlock->Block_List[j].Size==1) { lhs_rhs_done=true; @@ -1089,10 +1202,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } else lhs_rhs_done=false; - /*if (ModelBlock->Block_List[j].Size==1) - lhs_rhs_done=true; - else - lhs_rhs_done=false;*/ //The Temporary terms temporary_terms_type tt2; i=0; @@ -1124,7 +1233,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl // 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); + //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); if (!lhs_rhs_done) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; @@ -1343,6 +1452,8 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n"; #endif break; + default: + break; } prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type; @@ -1734,13 +1845,10 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size; int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j]; int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j]; - /*cout << " ! IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr+k1*block_triangular.ModelBlock->Block_List[num].Size << "), " << k1 << ")] = " << u << ";\n"; - cout << " ? IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr << "), " << k1 << ")] = " << u << ";\n";*/ SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&u), sizeof(u)); - //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n"; u_count_int++; } } @@ -1748,13 +1856,12 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b { int eqr1=j; int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods - +block_triangular.Model_Max_Lead); + +block_triangular.incidencematrix.Model_Max_Lead_Endo); int k1=0; SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); SaveCode.write(reinterpret_cast(&k1), sizeof(k1)); SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); - //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n"; u_count_int++; } for(j=0;jBlock_List[num].Size;j++) @@ -1775,7 +1882,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b { string filename; ofstream mStaticModelFile; - int i, k, prev_Simulation_Type; + int i, k, prev_Simulation_Type, ga_index = 1; bool /*printed = false,*/ skip_head, open_par=false; filename = static_basename + ".m"; mStaticModelFile.open(filename.c_str(), ios::out | ios::binary); @@ -1812,6 +1919,12 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; } mStaticModelFile << " ];\n"; + mStaticModelFile << " y_index_eq=["; + for(int ik=0;ikBlock_List[i].Size;ik++) + { + mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; + } + mStaticModelFile << " ];\n"; k=block_triangular.ModelBlock->Block_List[i].Simulation_Type; if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) @@ -1825,21 +1938,27 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b case EVALUATE_FORWARD_R: case EVALUATE_BACKWARD_R: if(!skip_head) - mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; + { + ga_index = 1; + mStaticModelFile << " [y, ga]=" << static_basename << "_" << i + 1 << "(y, x, 1);\n"; + } + else + ga_index++; mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n"; + mStaticModelFile << " g1(y_index_eq, y_index) = ga(" << ga_index << ", " << ga_index << ");\n"; break; case SOLVE_FORWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_SIMPLE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - mStaticModelFile << " [r, g1]=" << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1(y_index_eq, y_index)]=" << static_basename << "_" << i + 1 << "(y, x, 1);\n"; mStaticModelFile << " residual(y_index)=r;\n"; break; } prev_Simulation_Type=k; } - mStaticModelFile << " varargout{1}=residual;\n"; + mStaticModelFile << " varargout{1}=residual';\n"; mStaticModelFile << " varargout{2}=g1;\n"; mStaticModelFile << " return;\n"; mStaticModelFile << " end;\n"; @@ -1867,7 +1986,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b { mStaticModelFile << " end\n"; } - mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " y=" << static_basename << "_" << i + 1 << "(y, x, 0);\n"; } open_par=false; } @@ -1880,32 +1999,16 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b open_par=false; mStaticModelFile << " g1=0;\n"; mStaticModelFile << " r=0;\n"; - /*mStaticModelFile << " for it_=y_kmin+1:periods+y_kmin\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " Per_y_=it_*y_size;\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r/g1;\n"; - mStaticModelFile << " cvg=((r[0]*r[0])maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mStaticModelFile << " cvg=((r*r)maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n"; - mStaticModelFile << " [L, U] = LU(g1);\n"; - mStaticModelFile << " y(it_, y_index) = U\\(L\\b);\n"; - mStaticModelFile << " cvg=((r'*r)maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n"; + mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; mStaticModelFile << " max_res=max(abs(r));\n"; mStaticModelFile << " cvg=(max_resBlock_List[i].Size;ik++) { tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; } - //mDynamicModelFile << " ];\n"; if(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R) { if(i==block_triangular.ModelBlock->Size-1) @@ -2127,8 +2198,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_SIMPLE: mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n"; - mDynamicModelFile << " y_index = " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n"; - mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n"; + mDynamicModelFile << " [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n"; mDynamicModelFile << " residual(y_index_eq)=r;\n"; tmp_eq.str(""); tmp.str(""); @@ -2136,41 +2206,38 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string case SOLVE_FORWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1; - mDynamicModelFile << " y_index_c = y_index;\n"; - mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n"; - mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n"; - mDynamicModelFile << " end;\n"; - //tmp_i=variable_table.max_lag+variable_table.max_lead+1; - mDynamicModelFile << " ga = [];\n"; - mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n"; - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n"; - mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n"; - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n"; + mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n"; mDynamicModelFile << " residual(y_index_eq)=r;\n"; break; case SOLVE_TWO_BOUNDARIES_COMPLETE: case SOLVE_TWO_BOUNDARIES_SIMPLE: - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, y_size, 1);\n"; + int j; mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; - mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " y_index = ["; + for(j=0;jBlock_List[i].Size;ik++) + { + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr; + } + int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1; + for(j=0;jBlock_List[i].nb_exo;ik++) + mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr+symbol_table.endo_nbr*tmp_i; + mDynamicModelFile << " ];\n"; tmp.str(""); tmp_eq.str(""); - tmp_i=variable_table.max_lag+variable_table.max_lead+1; mDynamicModelFile << " ga = [];\n"; - mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n"; - mDynamicModelFile << " y_index_c = y_index;\n"; - tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1; - mDynamicModelFile << " for i=1:" << tmp_i-1 << ",\n"; - mDynamicModelFile << " y_index_c = [y_index_c (y_index+i*y_size)];\n"; - mDynamicModelFile << " end;\n"; - //mDynamicModelFile << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n"; - mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n"; + j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1) + + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1); + mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " << + block_triangular.ModelBlock->Block_List[i].Size*j << ");\n"; + tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1; + mDynamicModelFile << " [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_-" << variable_table.max_lag << ", 1, ga, g2, g3, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n"; if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead) - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga;\n"; + mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n"; else - mDynamicModelFile << " g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n"; + mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n"; mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n"; break; } @@ -2256,7 +2323,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " iter=0;\n"; mDynamicModelFile << " Per_y_=it_*y_size;\n"; mDynamicModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n"; + mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n"; mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mDynamicModelFile << " cvg=((r*r)maxit_),\n"; - mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n"; + mDynamicModelFile << " [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n"; mDynamicModelFile << " y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; mDynamicModelFile << " cvg=((r*r)Block_List[i].Size)) { - /*if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - if (!printed) - { - printed = true; - } - SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr); - Nb_SGE++; - mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/ if (open_par) mDynamicModelFile << " end\n"; open_par=false; @@ -2348,17 +2405,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " return;\n"; mDynamicModelFile << " end\n"; mDynamicModelFile << " end\n"; + /*cerr << "Not implemented block SOLVE_FORWARD_COMPLETE" << endl; exit(EXIT_FAILURE);*/ } else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size)) { - /*if (open_par) - mDynamicModelFile << " end\n"; - open_par=false; - SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr); - Nb_SGE++; - mDynamicModelFile << " Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/ if (open_par) mDynamicModelFile << " end\n"; open_par=false; @@ -2398,7 +2450,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string printed = true; } Nb_SGE++; - //cout << "new_SGE=" << new_SGE << "\n"; mDynamicModelFile << " cvg=0;\n"; mDynamicModelFile << " iter=0;\n"; @@ -2410,13 +2461,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string } mDynamicModelFile << " ];\n"; mDynamicModelFile << " Blck_size=" << block_triangular.ModelBlock->Block_List[i].Size << ";\n"; - /*mDynamicModelFile << " if(options_.simulation_method==2 | options_.simulation_method==3),\n"; - mDynamicModelFile << " [r, g1]= " << basename << "_static(y, x);\n"; - mDynamicModelFile << " [L1,U1] = lu(g1,1e-5);\n"; - mDynamicModelFile << " I = speye(periods);\n"; - mDynamicModelFile << " L1=kron(I,L1);\n"; - mDynamicModelFile << " U1=kron(I,U1);\n"; - mDynamicModelFile << " end;\n";*/ mDynamicModelFile << " y_kmin_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lag << ";\n"; mDynamicModelFile << " y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n"; mDynamicModelFile << " lambda=options_.slowc;\n"; @@ -2428,12 +2472,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size; mDynamicModelFile << " Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n"; mDynamicModelFile << " g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n"; - /*mDynamicModelFile << " cpath=path;\n"; - mDynamicModelFile << " addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n";*/ - mDynamicModelFile << " bicgstabh=@bicgstab;\n"; - //mDynamicModelFile << " path(cpath);\n"; mDynamicModelFile << sp << " reduced = 0;\n"; - //mDynamicModelFile << " functions(bicgstabh)\n"; if (!block_triangular.ModelBlock->Block_List[i].is_linear) { sp=" "; @@ -2443,7 +2482,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string { sp=""; } - mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n"; + mDynamicModelFile << sp << " [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, periods, 0, g1, g2, g3, y_kmin, Blck_size);\n"; mDynamicModelFile << sp << " g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n"; mDynamicModelFile << sp << " b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n"; mDynamicModelFile << sp << " if(~isreal(r))\n"; @@ -2525,7 +2564,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << sp << " flag1=1;\n"; mDynamicModelFile << sp << " while(flag1>0)\n"; mDynamicModelFile << sp << " [L1, U1]=luinc(g1a,luinc_tol);\n"; - mDynamicModelFile << sp << " [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; + mDynamicModelFile << sp << " [za,flag1] = bicgstab(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n"; mDynamicModelFile << sp << " if (flag1>0 | reduced)\n"; mDynamicModelFile << sp << " if(flag1==1)\n"; mDynamicModelFile << sp << " disp(['No convergence inside BICGSTAB after ' num2str(periods*" << block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n"; @@ -2791,11 +2830,12 @@ ModelTree::writeOutput(ostream &output) const */ output << "M_.lead_lag_incidence = ["; // Loop on endogenous variables + int lag = 0; for(int endoID = 0; endoID < symbol_table.endo_nbr; endoID++) { output << "\n\t"; // Loop on periods - for(int lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) + for(lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++) { // Print variableID if exists with current period, otherwise print 0 try @@ -2819,7 +2859,7 @@ ModelTree::writeOutput(ostream &output) const bool skip_the_head; int k=0; int count_lead_lag_incidence = 0; - int max_lead, max_lag; + int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo; for(int j = 0;j < block_triangular.ModelBlock->Size;j++) { //For a block composed of a single equation determines wether we have to evaluate or to solve the equation @@ -2835,9 +2875,17 @@ ModelTree::writeOutput(ostream &output) const k++; count_lead_lag_incidence = 0; int Block_size=block_triangular.ModelBlock->Block_List[j].Size; + //int Block_exo_nbr=block_triangular.ModelBlock->Block_List[j].exo_nbr; max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ; max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead; + max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ; + max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo; + max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ; + max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo; bool evaluate=false; + vector exogenous;//, tmp_exogenous(block_triangular.exo_nbr,-1); + vector::iterator it_exogenous; + exogenous.clear(); ostringstream tmp_s, tmp_s_eq; tmp_s.str(""); tmp_s_eq.str(""); @@ -2846,6 +2894,13 @@ ModelTree::writeOutput(ostream &output) const tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1; tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1; } + for(int i=0;iBlock_List[j].nb_exo;i++) + { + int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i]; + for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/; + 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 @@ -2864,78 +2919,74 @@ ModelTree::writeOutput(ostream &output) const max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ; if(max_leadBlock_List[j+Block_size].Max_Lead) max_lead=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead; - //cout << "block_triangular.ModelBlock->Block_List[" << j+Block_size << "].Size=" << block_triangular.ModelBlock->Block_List[j+Block_size].Size << "\n"; + if(max_lag_endo 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_endoBlock_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_List[j+Block_size].Max_Lag_Exo ) + max_lag_exo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo ; + if(max_lead_exoBlock_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;iBlock_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;iBlock_List[j+Block_size].nb_exo;i++) + { + int ii=block_triangular.ModelBlock->Block_List[j+Block_size].Exogenous[i]; + //for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) cout << "*it_exogenous=" << *it_exogenous << "\n"; + if(it_exogenous==exogenous.end()) + exogenous.push_back(ii); + } Block_size+=block_triangular.ModelBlock->Block_List[j+Block_size].Size; } - //cout << "i=" << i << " max_lag=" << max_lag << " max_lead=" << max_lead << "\n"; } } output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n"; - //output << "M_.block_structure.block(" << k << ").size = " << block_triangular.ModelBlock->Block_List[j].Size << ";\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_endo_lag = " << max_lag << ";\n"; - output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n"; + output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n"; output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n"; output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n"; output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n"; + output << "M_.block_structure.block(" << k << ").exogenous = ["; + int i=0; + for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++) + if(*it_exogenous>=0) + { + output << " " << *it_exogenous+1; + i++; + } + output << "];\n"; + output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n"; + + output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n"; + tmp_s.str(""); - cout << "begining of lead_lag_incidence\n"; bool done_IM=false; if(!evaluate) { output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n"; - for(int l=-max_lag;lBlock_List[j].Size;l_var++) + tmp_IM=block_triangular.incidencematrix.bGet_IM(l, eEndogenous); + if(tmp_IM) { - for(int l_equ=0;l_equBlock_List[j].Size;l_equ++) - if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) - { - count_lead_lag_incidence++; - if(tmp_s.str().length()) - tmp_s << " "; - tmp_s << count_lead_lag_incidence; - done_IM=true; - break; - } - if(!done_IM) - tmp_s << " 0"; - done_IM=false; - } - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; - tmp_s.str(""); - } - } - else - { - output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; - for(int l=-max_lag;lBlock_List[ii].Size;l_var++) + for(int l_var=0;l_varBlock_List[j].Size;l_var++) { - for(int l_equ=0;l_equBlock_List[ii].Size;l_equ++) - if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + for(int l_equ=0;l_equBlock_List[j].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]]) { - //if(not_increm && l==-max_lag) - count_lead_lag_incidence++; - not_increm=false; + count_lead_lag_incidence++; if(tmp_s.str().length()) tmp_s << " "; - //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; tmp_s << count_lead_lag_incidence; done_IM=true; break; @@ -2944,10 +2995,53 @@ ModelTree::writeOutput(ostream &output) const tmp_s << " 0"; done_IM=false; } - ii++; + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n"; + tmp_s.str(""); + } + } + } + else + { + bool done_some_where; + output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n"; + for(int l=-max_lag_endo;lBlock_List[ii].Size;l_var++) + { + for(int l_equ=0;l_equBlock_List[ii].Size;l_equ++) + if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]]) + { + //if(not_increm && l==-max_lag) + count_lead_lag_incidence++; + not_increm=false; + if(tmp_s.str().length()) + tmp_s << " "; + //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size; + tmp_s << count_lead_lag_incidence; + done_IM=true; + break; + } + if(!done_IM) + tmp_s << " 0"; + else + done_some_where = true; + done_IM=false; + } + ii++; + } + //if(done_some_where) + output << tmp_s.str() << "\n"; + tmp_s.str(""); } - output << tmp_s.str() << "\n"; - tmp_s.str(""); } output << "];\n"; } @@ -2955,13 +3049,13 @@ ModelTree::writeOutput(ostream &output) const prev_Simulation_Type=block_triangular.ModelBlock->Block_List[j].Simulation_Type; } - for(int j=-block_triangular.Model_Max_Lag;jfirst.second) == eEndogenous) { NodeID Id = it->second; - double val; + double val = 0; try { val = Id->eval(eval_context); @@ -3046,7 +3140,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_ int k1=variable_table.getLag(it->first.second); if (a_variable_lag!=k1) { - IM=block_triangular.bGet_IM(k1); + IM=block_triangular.incidencematrix.bGet_IM(k1, eEndogenous); a_variable_lag=k1; } if (k1==0) @@ -3058,7 +3152,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_ { if(block_triangular.bt_verbose) cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr << ")\n"; - block_triangular.unfill_IM(eq, var, k1); + block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous); i++; } } @@ -3159,6 +3253,8 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term if (mode == eSparseDLLMode || mode == eSparseMode) { + block_triangular.incidencematrix.init_incidence_matrix(); + BuildIncidenceMatrix(); jacob_map j_m; evaluateJacobian(eval_context, &j_m); @@ -3166,7 +3262,7 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term if (block_triangular.bt_verbose) { cout << "The gross incidence matrix \n"; - block_triangular.Print_IM( symbol_table.endo_nbr); + block_triangular.incidencematrix.Print_IM(eEndogenous); } block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m); BlockLinear(block_triangular.ModelBlock); @@ -3208,7 +3304,8 @@ ModelTree::writeDynamicFile(const string &basename) const case eSparseMode: writeSparseDynamicMFile(basename + "_dynamic", basename, mode); block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.Free_IM(block_triangular.First_IM); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); break; case eDLLMode: writeDynamicCFile(basename + "_dynamic"); @@ -3216,7 +3313,8 @@ ModelTree::writeDynamicFile(const string &basename) const case eSparseDLLMode: writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL); block_triangular.Free_Block(block_triangular.ModelBlock); - block_triangular.Free_IM(block_triangular.First_IM); + block_triangular.incidencematrix.Free_IM(); + //block_triangular.Free_IM_X(block_triangular.First_IM_X); break; } } diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index 1516a026a..c1109e1e7 100644 --- a/preprocessor/ParsingDriver.cc +++ b/preprocessor/ParsingDriver.cc @@ -169,7 +169,7 @@ ParsingDriver::add_model_variable(string *name, string *olag) if (type == eUnknownFunction) error("Symbol " + *name + " is a function name unknown to Dynare. It cannot be used inside model."); - if (type == eExogenous && lag != 0) + if (type == eExogenous && lag != 0 && (model_tree->mode != eSparseDLLMode && model_tree->mode != eSparseMode)) warning("Exogenous variable " + *name + " has lead/lag " + *olag); if (type == eModelLocalVariable && lag != 0) @@ -179,9 +179,13 @@ ParsingDriver::add_model_variable(string *name, string *olag) NodeID id = model_tree->AddVariable(*name, lag); - if ((type == eEndogenous) && (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode)) - model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); - + /*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode) + { + if (type == eEndogenous) + model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); + if (type == eExogenous) + model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); + }*/ delete name; delete olag; return id; @@ -360,14 +364,16 @@ void ParsingDriver::sparse_dll() { model_tree->mode = eSparseDLLMode; - model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/ } void ParsingDriver::sparse() { model_tree->mode = eSparseMode; - model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr); + model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/ } void diff --git a/preprocessor/include/BlockTriangular.hh b/preprocessor/include/BlockTriangular.hh index 894b8ccbd..fd9e7898c 100644 --- a/preprocessor/include/BlockTriangular.hh +++ b/preprocessor/include/BlockTriangular.hh @@ -36,41 +36,61 @@ struct List_IM bool* IM; }; + +//! create and manage the incidence matrix +class IncidenceMatrix //: SymbolTable +{ + //friend class BlockTriangular; +public: +const SymbolTable &symbol_table; + IncidenceMatrix(const SymbolTable &symbol_table_arg); + List_IM* Build_IM(int lead_lag, SymbolType type); + List_IM* Get_IM(int lead_lag, SymbolType type) const; + bool* bGet_IM(int lead_lag, SymbolType type) const; + void fill_IM(int equation, int variable_endo, int lead_lag, SymbolType type); + void unfill_IM(int equation, int variable_endo, int lead_lag, SymbolType type); + void init_incidence_matrix(); + void Free_IM() const; + List_IM* Get_First(SymbolType type) const; + void Print_IM(SymbolType type) const; + void Print_SIM(bool* IM, SymbolType type) const; + + void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const; +private: + List_IM *First_IM, *Last_IM, *First_IM_X, *Last_IM_X ; +public: + int Model_Max_Lead, Model_Max_Lag; + int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo; +}; + + //! Matrix of doubles for representing jacobian typedef map,double> jacob_map; //! Create the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition class BlockTriangular { + //friend class IncidenceMatrix; public: - BlockTriangular(const SymbolTable &symbol_table_arg); const SymbolTable &symbol_table; + BlockTriangular(const SymbolTable &symbol_table_arg); + //BlockTriangular(const IncidenceMatrix &incidence_matrix_arg); + //const SymbolTable &symbol_table; Blocks blocks; Normalization normalization; - List_IM* Build_IM(int lead_lag); - List_IM* Get_IM(int lead_lag); - bool* bGet_IM(int lead_lag) const; - void fill_IM(int equation, int variable_endo, int lead_lag); - void unfill_IM(int equation, int variable_endo, int lead_lag); - void init_incidence_matrix(int nb_endo); - void Print_IM(int n) const; - void Free_IM(List_IM* First_IM) const; - void Print_SIM(bool* IM, int n) const; + IncidenceMatrix incidencematrix; void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m); bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m); void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0); - void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n); void Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock); void Free_Block(Model_Block* ModelBlock) const; - List_IM *First_IM ; - List_IM *Last_IM ; simple *Index_Equ_IM; simple *Index_Var_IM; int prologue, epilogue; - int Model_Max_Lead, Model_Max_Lag, periods; bool bt_verbose; - int endo_nbr; + //int endo_nbr, exo_nbr; Model_Block* ModelBlock; + int periods; inline static std::string BlockType0(int type) { switch (type) @@ -98,14 +118,14 @@ public: { case EVALUATE_FORWARD: case EVALUATE_FORWARD_R: - return ("EVALUATE FORWARD "); + return ("EVALUATE FORWARD "); break; case EVALUATE_BACKWARD: case EVALUATE_BACKWARD_R: return ("EVALUATE BACKWARD "); break; case SOLVE_FORWARD_SIMPLE: - return ("SOLVE FORWARD SIMPLE "); + return ("SOLVE FORWARD SIMPLE "); break; case SOLVE_BACKWARD_SIMPLE: return ("SOLVE BACKWARD SIMPLE "); @@ -114,7 +134,7 @@ public: return ("SOLVE TWO BOUNDARIES SIMPLE "); break; case SOLVE_FORWARD_COMPLETE: - return ("SOLVE FORWARD COMPLETE "); + return ("SOLVE FORWARD COMPLETE "); break; case SOLVE_BACKWARD_COMPLETE: return ("SOLVE BACKWARD COMPLETE "); diff --git a/preprocessor/include/ExprNode.hh b/preprocessor/include/ExprNode.hh index 55a2f9e12..d3d24e351 100644 --- a/preprocessor/include/ExprNode.hh +++ b/preprocessor/include/ExprNode.hh @@ -143,6 +143,7 @@ public: /*! Endogenous are stored as integer pairs of the form (symb_id, lag) They are added to the set given in argument */ virtual void collectEndogenous(set > &result) const = 0; + virtual void collectExogenous(set > &result) const = 0; virtual void computeTemporaryTerms(map &reference_count, temporary_terms_type &temporary_terms, map &first_occurence, @@ -179,6 +180,7 @@ public: NumConstNode(DataTree &datatree_arg, int id_arg); virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -198,6 +200,7 @@ public: VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg); virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -223,6 +226,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -250,6 +254,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -282,6 +287,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException); virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -307,6 +313,7 @@ public: Model_Block *ModelBlock, map_idx_type &map_idx) const; virtual void collectEndogenous(set > &result) const; + virtual void collectExogenous(set > &result) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; }; @@ -314,21 +321,22 @@ public: //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model struct IM_compact { - int size, u_init, u_finish, nb_endo; - int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index; + int size, u_init, u_finish, nb_endo, size_exo; + int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index; }; //! One block of the model struct Block { - int Size, Sized; + int Size, Sized, nb_exo, nb_exo_det; BlockType Type; BlockSimulationType Simulation_Type; int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo; + int Max_Lag_Endo, Max_Lead_Endo; + int Max_Lag_Exo, Max_Lead_Exo; bool is_linear; int *Equation, *Own_Derivative; - int *Variable, *Variable_Sorted, *dVariable; - int *variable_dyn_index, *variable_dyn_leadlag; + int *Variable, *Exogenous; temporary_terms_type *Temporary_terms; IM_compact *IM_lead_lag; int Code_Start, Code_Length; @@ -339,7 +347,7 @@ struct Model_Block { int Size, Periods; Block* Block_List; - int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block; + //int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block; }; #endif diff --git a/preprocessor/include/ModelTree.hh b/preprocessor/include/ModelTree.hh index 346455999..058bb9175 100644 --- a/preprocessor/include/ModelTree.hh +++ b/preprocessor/include/ModelTree.hh @@ -26,6 +26,7 @@ using namespace std; #include #include #include +#include #include "SymbolTable.hh" #include "NumericalConstants.hh" @@ -82,12 +83,14 @@ private: //! Computes derivatives of ModelTree void derive(int order); //! Write derivative of an equation w.r. to a variable - void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; + void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) const; //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const; //! Computes temporary terms void computeTemporaryTerms(int order); void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock); + //! Build The incidence matrix form the modeltree + void BuildIncidenceMatrix(); //! Writes temporary terms void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const; //! Writes model local variables diff --git a/tests/ferhat/fs2000.mod b/tests/ferhat/fs2000.mod index c223ba774..673598ac8 100644 --- a/tests/ferhat/fs2000.mod +++ b/tests/ferhat/fs2000.mod @@ -48,7 +48,7 @@ e = exp(e_a); y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); gy_obs = dA*y/y(-1); gp_obs = (P/P(-1))*m(-1)/dA; -vv = 0.2*ww+0.5*vv(-1)+1; +vv = 0.2*ww+0.5*vv(-1)+1+c(-1)+e_a; ww = 0.1*vv+0.5*ww(-1)+2; /* A lt= 0.5*vv-0.2*ww = 1 @@ -84,11 +84,11 @@ vv = 0; ww = 0; end; -/*shocks; +shocks; var e_a; stderr 0.014; var e_m; stderr 0.005; end; -*/ + options_.solve_tolf=1e-10; options_.maxit_=100; steady; @@ -102,6 +102,8 @@ end; simul(periods=200, method=lu); +stoch_simul(periods=200,order=1); + rplot y; rplot k; rplot c; diff --git a/tests/ferhat/ireland.mod b/tests/ferhat/ireland.mod index a7d0604a4..c51cba745 100644 --- a/tests/ferhat/ireland.mod +++ b/tests/ferhat/ireland.mod @@ -89,7 +89,7 @@ end; options_.maxit_=20; model_info; -simul(periods=200, method=/*LU*/GMRES/*bicgstab*/); +simul(periods=2000, method=/*LU*/GMRES/*bicgstab*/); rplot y; rplot k; diff --git a/tests/ferhat/ls2003.mod b/tests/ferhat/ls2003.mod index 21e2cdf62..8a8cb8238 100644 --- a/tests/ferhat/ls2003.mod +++ b/tests/ferhat/ls2003.mod @@ -1,4 +1,4 @@ -var y y_s R pie dq pie_s de A y_obs pie_obs R_obs; +var y y_s R pie dq pie_s de A y_obs pie_obs R_obs vv ww; varexo e_R e_q e_ys e_pies e_A; parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies; @@ -31,9 +31,22 @@ A = rho_A*A(-1)+e_A; y_obs = y-y(-1)+A; pie_obs = 4*pie; R_obs = 4*R; +vv = 0.2*ww+0.5*vv(-1)+1; +ww = 0.1*vv+0.5*ww(-1)+2; +/* A lt= + 0.5*vv-0.2*ww = 1 +-0.1*vv+0.5*ww = 2 +[ 0.5 -0.2][vv] [1] + = +[-0.1 0.5][ww] [2] +det = 0.25-0.02 = 0.23 +[vv] [0.5 0.2] [1] [0.9] [3.91304] + = 1/0.23* = 1/0.23* = +[ww] [0.1 0.5] [2] [1.1] [4.7826] +*/ end; -/*shocks; +shocks; var e_R = 1.25^2; var e_q = 2.5^2; var e_A = 1.89; @@ -41,7 +54,7 @@ var e_ys = 1.89; var e_pies = 1.89; end; -varobs y_obs R_obs pie_obs dq de; +/*varobs y_obs R_obs pie_obs dq de; estimated_params; psi1 , gamma_pdf,1.5,0.5; @@ -63,12 +76,13 @@ stderr e_ys,inv_gamma_pdf,1.2533,0.6551; stderr e_pies,inv_gamma_pdf,1.88,0.9827; end; -estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=10,prefilter=1,mh_jscale=0.5,mh_replic=0); - +estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=10,prefilter=1,mh_jscale=0.5,mh_replic=0,nograph); */ + +options_.maxit_=100; steady; -//model_info; +model_info; check; shocks; @@ -77,6 +91,9 @@ periods 1; values 0.5; end; -//simul(periods=200,method=bicgstab); -//rplot A; -//rplot pie; +simul(periods=200,method=bicgstab); +rplot A; +rplot pie; + +stoch_simul(periods=200,order=1); + diff --git a/tests/ferhat/ramst.mod b/tests/ferhat/ramst.mod index 8285d2810..8e64dbf04 100644 --- a/tests/ferhat/ramst.mod +++ b/tests/ferhat/ramst.mod @@ -32,7 +32,7 @@ periods 1; values 1.002; end; -simul(periods=200, METHOD=GmRes); +simul(periods=200, METHOD=LU); rplot c; rplot k;