diff --git a/matlab/check.m b/matlab/check.m
index 457eb9d1d88584e6e29ab9ef250b1277f7f2c349..15642a22cf6525fe84fcbd1e10cbf96786367d16 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -59,7 +59,7 @@ global it_
if(options_.model_mode == 1)
nyf = dr.nyf;
else
- nyf = nnz(dr.kstate(:,2)>M_.maximum_lag+1);
+ nyf = nnz(dr.kstate(:,2)>M_.maximum_endo_lag+1);
end;
[m_lambda,i]=sort(abs(eigenvalues_));
n_explod = nnz(abs(eigenvalues_) > options_.qz_criterium);
diff --git a/matlab/dr11_sparse.m b/matlab/dr11_sparse.m
index 7fddf48e7a37bfe70e85760ff262d12ad3e9e65a..1a2f7cc0528b9420c1a76a28c8e021e689ae2124 100644
--- a/matlab/dr11_sparse.m
+++ b/matlab/dr11_sparse.m
@@ -23,34 +23,18 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
kstate = dr.kstate;
kad = dr.kad;
kae = dr.kae;
- %kstate
- %kad
- %kae
nstatic = dr.nstatic;
nfwrd = dr.nfwrd;
npred = dr.npred;
nboth = dr.nboth;
- %nstatic
- %nfwrd
- %npred
- %nboth
order_var = dr.order_var;
- %order_var
nd = size(kstate,1);
- %nd
nz = nnz(M_.lead_lag_incidence);
- %nz
sdyn = M_.endo_nbr - nstatic;
- %sdyn
-% M_.lead_lag_incidence
k0 = M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var);
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
-% size(jacobia_)
-% k0
- %k1
b = jacobia_(:,k0);
- %full(b)
if M_.maximum_endo_lead == 0; % backward models
a = jacobia_(:,nonzeros(k1'));
@@ -75,7 +59,6 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
end
return;
end
-
%forward--looking models
if nstatic > 0
[Q,R] = qr(b(:,1:nstatic));
@@ -83,14 +66,8 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
else
aa = jacobia_;
end
-% full(aa)
a = aa(:,nonzeros(k1'));
b = aa(:,k0);
- %M_.lead_lag_incidence
- %k0
- %k1
- %a
- %b
b10 = b(1:nstatic,1:nstatic);
b11 = b(1:nstatic,nstatic+1:end);
b2 = b(nstatic+1:end,nstatic+1:end);
@@ -100,9 +77,9 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
end
% buildind D and E
+ %nd
d = zeros(nd,nd) ;
e = d ;
-
k = find(kstate(:,2) >= M_.maximum_endo_lag+2 & kstate(:,3));
d(1:sdyn,k) = a(nstatic+1:end,kstate(k,3)) ;
k1 = find(kstate(:,2) == M_.maximum_endo_lag+2);
@@ -124,10 +101,6 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
[ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options_.qz_criterium);
- %ss
- %tt
- %sdim
- %fprintf('%20.16f\n',dr.eigval)
if info1
info(1) = 2;
@@ -138,14 +111,13 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
nba = nd-sdim;
nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
- %disp(['task' num2str(task)]);
+
if task == 1
dr.rank = rank(w(1:nyf,nd-nyf+1:end));
% Under Octave, eig(A,B) doesn't exist, and
% lambda = qz(A,B) won't return infinite eigenvalues
if ~exist('OCTAVE_VERSION')
dr.eigval = eig(e,d);
-% dr.eigval
end
return
end
diff --git a/matlab/dr1_sparse.m b/matlab/dr1_sparse.m
index 064d307ceb2f5cbbc04c028fe3865577b540c7fa..acef2657e99b6cf05a099c65259ac7dcf96411a3 100644
--- a/matlab/dr1_sparse.m
+++ b/matlab/dr1_sparse.m
@@ -214,11 +214,12 @@ function [dr,info,M_,options_,oo_] = dr1_sparse(dr,task,M_,options_,oo_)
dr.nyf = 0;
dr.rank = 0;
for i=1:length(M_.block_structure.block)
- %disp(['block ' num2str(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));
- jcb_=jacobia_(M_.block_structure.block(i).equation,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)));
- jcb_=jcb_(:,find(any(jcb_,1)));
+ 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;
+ jcb_=jacobia_(row_selector,col_selector);
+ jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ;
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;
diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index a400646cdd77e1b960d405df9b4aef077dfde04b..f17104b7477ad53dd96df63b255e4d9797c4cbd5 100755
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 7a110a4ca84d2c049ecc724a398244b357139bd0..c6bef42d5faa37f0f5e1d691bb5684355186a88f 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -32,7 +32,7 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
DataTree(symbol_table_arg, num_constants_arg),
mode(eStandardMode),
- cutoff(1e-12),
+ cutoff(1e-15),
markowitz(0.7),
new_SGE(true),
computeJacobian(false),
@@ -543,7 +543,9 @@ 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(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 << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
<< "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
@@ -2124,9 +2126,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
break;
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_SIMPLE:
- mDynamicModelFile << " y_index=" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n";
+ 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 << " residual(y_index_eq)=r;\n";
+ tmp_eq.str("");
+ tmp.str("");
break;
case SOLVE_FORWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
@@ -2140,6 +2145,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
//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 << " residual(y_index_eq)=r;\n";
@@ -2159,6 +2165,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
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";
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";
@@ -2877,6 +2884,8 @@ ModelTree::writeOutput(ostream &output) const
output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n";
tmp_s.str("");
+ cout << "begining of lead_lag_incidence\n";
+
bool done_IM=false;
if(!evaluate)
{
diff --git a/tests/ferhat/fs2000.mod b/tests/ferhat/fs2000.mod
index 645db083366ca2ed2067e67f3afae75d4c1326d3..c223ba774233ef5fa7075c91a3b54771a16ecab4 100644
--- a/tests/ferhat/fs2000.mod
+++ b/tests/ferhat/fs2000.mod
@@ -18,7 +18,7 @@
//
// Michel Juillard, February 2004
-var m P c e W R k d n l gy_obs gp_obs y dA;
+var m P c e W R k d n l gy_obs gp_obs y dA vv ww;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
@@ -48,6 +48,19 @@ 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;
+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;
initval;
@@ -67,6 +80,8 @@ gp_obs = exp(-gam);
dA = exp(gam);
e_a=0;
e_m=0;
+vv = 0;
+ww = 0;
end;
/*shocks;
diff --git a/tests/ferhat/ls2003.mod b/tests/ferhat/ls2003.mod
index c8339f5a87ba60b3439f718d6301d4f02beed0ec..21e2cdf626c220658248caf560ce8dd1b1850215 100644
--- a/tests/ferhat/ls2003.mod
+++ b/tests/ferhat/ls2003.mod
@@ -17,7 +17,7 @@ rho_ys = 0.9;
rho_pies = 0.7;
-//model(sparse_dll,gcc_compiler,cutoff=1e-17);
+//model(sparse_dll,cutoff=1e-17);
model(sparse);
//model;
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
@@ -68,7 +68,7 @@ estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=10,prefilter=1,mh_js
*/
steady;
-model_info;
+//model_info;
check;
shocks;
@@ -77,6 +77,6 @@ periods 1;
values 0.5;
end;
-simul(periods=200,method=bicgstab);
-rplot A;
-rplot pie;
+//simul(periods=200,method=bicgstab);
+//rplot A;
+//rplot pie;