Commit 7c6ed357 authored by ferhat's avatar ferhat
Browse files

Bugs correction in check command with sparse option in model command

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2196 ac1d8469-bf42-47a9-8791-bf33cf982152
parent dcc37a68
......@@ -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);
......
......@@ -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
......
......@@ -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;
......
No preview for this file type
......@@ -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)
{
......
......@@ -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;
......
......@@ -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;
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment