Commit e85af5dc authored by adjemian's avatar adjemian
Browse files

Added a dynare configuration routine and bug corrections.

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1826 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 4f50ae6e
......@@ -11,7 +11,7 @@ function [dr,info]=dr1(dr,task)
% OUTPUTS
% dr: structure of decision rules for stochastic simulations
% info = 1: the model doesn't define current variables uniquely
% info = 2: problem in mjdgges.dll info(2) contains error code
% info = 2: problem in mjdgges.dll info(2) contains error code.
% info = 3: BK order condition not satisfied info(2) contains "distance"
% absence of stable trajectory
% info = 4: BK order condition not satisfied info(2) contains "distance"
......@@ -266,8 +266,6 @@ if ~isempty(kad)
end
end
use_qzdiv = 0;
% $$$ if exist('ordqz')
% $$$ info1 = 0;
% $$$ try
......@@ -284,28 +282,23 @@ use_qzdiv = 0;
% $$$ end
% $$$ nba = nd-sdim;
% $$$ elseif exist('mjdgges')
if exist('mjdgges')
[ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options_.qz_criterium);
if info1
if exist('mjdgges')==2
use_qzdiv = 1;
else
use_qzdiv = 0;
end
[ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options_.qz_criterium);
if info1
info(1) = 2;
info(2) = info1;
return
end
nba = nd-sdim;
else
% using Chris Sim's routines
use_qzdiv = 1;
[ss,tt,qq,w] = qz(e,d);
[tt,ss,qq,w] = qzdiv(options_.qz_criterium,tt,ss,qq,w);
ss1=diag(ss);
tt1=diag(tt);
warning_state = warning;
warning off;
dr.eigval = ss1./tt1 ;
warning warning_state;
nba = nnz(abs(dr.eigval) > options_.qz_criterium);
end
nba = nd-sdim;
nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
if task == 1
......@@ -387,7 +380,7 @@ if options_.loglinear == 1
dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
end
%necessary when using Sims' routines
%% Necessary when using Sims' routines for QZ
if use_qzdiv
gx = real(gx);
hx = real(hx);
......@@ -467,23 +460,7 @@ zx=[zx; zeros(M_.exo_nbr,np);zeros(M_.exo_det_nbr,np)];
zu=[zu; eye(M_.exo_nbr);zeros(M_.exo_det_nbr,M_.exo_nbr)];
[nrzx,nczx] = size(zx);
if ~exist('sparse_hessian_times_B_kronecker_C')
if nrzx*nrzx*nczx*nczx > 1e7
rhs = zeros(M_.endo_nbr,nczx*nczx);
k1 = 1;
for i1 = 1:nczx
for i2 = 1:nczx
rhs(:,k1) = hessian*kron(zx(:,i1),zx(:,i2));
k1 = k1 + 1;
end
end
else
rhs = hessian*kron(zx,zx);
end
else
rhs = sparse_hessian_times_B_kronecker_C(hessian,zx);
end
rhs = -rhs;
rhs = -sparse_hessian_times_B_kronecker_C(hessian,zx);
%lhs
n = M_.endo_nbr+sum(kstate(:,2) > M_.maximum_endo_lag+1 & kstate(:,2) < M_.maximum_endo_lag+M_.maximum_endo_lead+1);
......@@ -530,13 +507,8 @@ A(1:M_.endo_nbr,nstatic+1:nstatic+npred)=...
C = hx;
D = [rhs; zeros(n-M_.endo_nbr,size(rhs,2))];
if exist('gensylv')
dr.ghxx = gensylv(2,A,B,C,D);
else
C = kron(hx,hx);
x0 = sylvester3(A,B,C,D);
dr.ghxx = sylvester3a(x0,A,B,C,D);
end
dr.ghxx = gensylv(2,A,B,C,D);
%ghxu
%rhs
......@@ -544,44 +516,16 @@ hu = dr.ghu(nstatic+1:nstatic+npred,:);
%kk = reshape([1:np*np],np,np);
%kk = kk(1:npred,1:npred);
%rhs = -hessian*kron(zx,zu)-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hx(1:npred,:),hu(1:npred,:));
if ~exist('sparse_hessian_times_B_kronecker_C')
if nrzx*nrzx*nczx*M_.exo_nbr > 1e7
rhs = zeros(M_.endo_nbr,nczx*M_.exo_nbr);
k1 = 1;
for i1 = 1:nczx
for i2 = 1:M_.exo_nbr
rhs(:,k1) = hessian*kron(zx(:,i1),zu(:,i2));
k1 = k1 + 1;
end
end
else
rhs = hessian*kron(zx,zu);
end
else
rhs = sparse_hessian_times_B_kronecker_C(hessian,zx,zu);
end
rhs = sparse_hessian_times_B_kronecker_C(hessian,zx,zu);
nyf1 = sum(kstate(:,2) == M_.maximum_endo_lag+2);
hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
%B1 = [B(1:M_.endo_nbr,:);zeros(size(A,1)-M_.endo_nbr,size(B,2))];
[nrhx,nchx] = size(hx);
[nrhu1,nchu1] = size(hu1);
if ~exist('A_times_B_kronecker_C')
if nrhx*nrhu1*nchx*nchu1 > 1e7
B1 = zeros(size(dr.ghxx,1),nchx*nchu1);
k1 = 1;
for i1 = 1:nchx
for i2 = 1:nchu1
B1(:,k1) = dr.ghxx*kron(hx(:,i1),hu1(:,i2));
k1 = k1 + 1;
end
end
B1 = B*B1;
else
B1 = B*dr.ghxx*kron(hx,hu1);
end
else
B1 = B*A_times_B_kronecker_C(dr.ghxx,hx,hu1);
end
B1 = B*A_times_B_kronecker_C(dr.ghxx,hx,hu1);
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
......@@ -592,39 +536,11 @@ dr.ghxu = A\rhs;
%rhs
kk = reshape([1:np*np],np,np);
kk = kk(1:npred,1:npred);
if ~exist('sparse_hessian_times_B_kronecker_C')
if nrzx*nrzx*M_.exo_nbr*M_.exo_nbr > 1e7
rhs = zeros(M_.endo_nbr,M_.exo_nbr*M_.exo_nbr);
k1 = 1;
for i1 = 1:M_.exo_nbr
for i2 = 1:M_.exo_nbr
rhs(:,k1) = hessian*kron(zu(:,i1),zu(:,i2));
k1 = k1 + 1;
end
end
else
rhs = hessian*kron(zu,zu);
end
else
rhs = sparse_hessian_times_B_kronecker_C(hessian,zu);
end
if ~exist('A_times_B_kronecker_C')
if nrhu1*nrhu1*nchu1*nchu1 > 1e7
B1 = zeros(size(dr.ghxx,1),nchu1*nchu1);
k1 = 1;
for i1 = 1:nchu1
for i2 = 1:nchu1
B1(:,k1) = dr.ghxx*kron(hu1(:,i1),hu1(:,i2));
k1 = k1 + 1;
end
end
B1 = B*B1;
else
B1 = B*dr.ghxx*kron(hu1,hu1);
end
else
B1 = A_times_B_kronecker_C(B*dr.ghxx,hu1);
end
rhs = sparse_hessian_times_B_kronecker_C(hessian,zu);
B1 = A_times_B_kronecker_C(B*dr.ghxx,hu1);
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
%lhs
......@@ -668,27 +584,9 @@ for i=1:M_.maximum_endo_lead
[junk,k3a,k3] = ...
find(M_.lead_lag_incidenceordered(M_.maximum_endo_lag+j+1,:));
nk3a = length(k3a);
if ~exist('sparse_hessian_times_B_kronecker_C')
if nk3a*nk3a*M_.exo_nbr*M_.exo_nbr > 1e7
B1 = zeros(M_.endo_nbr,M_.exo_nbr*M_.exo_nbr);
k1 = 1;
Hesse = hessian(:,kh(k3,k3));
guk3a = gu(k3a,:);
for i1 = 1:M_.exo_nbr
for i2 = 1:M_.exo_nbr
B1(:,k1) = Hesse*kron(guk3a(:,i1),guk3a(:,i2));
k1 = k1 + 1;
end
end
else
B1 = hessian(:,kh(k3,k3))*kron(gu(k3a,:),gu(k3a,:));
end
else
B1 = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:));
end
B1 = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:));
RHS = RHS + jacobia_(:,k2)*guu(k2a,:)+B1;
end
% LHS
[junk,k2a,k2] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+i+1,order_var));
LHS = LHS + jacobia_(:,k2)*(E(k2a,:)+[O1(k2a,:) dr.ghx(k2a,:)*H O2(k2a,:)]);
......@@ -700,28 +598,8 @@ for i=1:M_.maximum_endo_lead
kk = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
gu = dr.ghx*Gu;
[nrGu,ncGu] = size(Gu);
if ~exist('A_times_B_kronecker_C')
if nrGu*nrGu*ncGu*ncGu > 1e7
G1 = zeros(M_.endo_nbr,ncGu*ncGu);
G2 = zeros(size(hxx,1),ncGu*ncGu);
k1 = 1;
for i1 = 1:nchx
for i2 = 1:nchu1
GuGu = kron(Gu(:,i1),Gu(:,i2));
G1(:,k1) = dr.ghxx*GuGu;
G2(:,k1) = hxx*GuGu;
k1 = k1 + 1;
end
end
else
GuGu = kron(Gu,Gu);
G1 = dr.ghxx*GuGu;
G2 = hxx*GuGu;
end
else
G1 = A_times_B_kronecker_C(dr.ghxx,Gu);
G2 = A_times_B_kronecker_C(hxx,Gu);
end
G1 = A_times_B_kronecker_C(dr.ghxx,Gu);
G2 = A_times_B_kronecker_C(hxx,Gu);
guu = dr.ghx*Guu+G1;
Gu = hx*Gu;
Guu = hx*Guu;
......
function dynare(fname, varargin)
% function dynare(fname, varargin)
% This command runs dynare with specified model file in argument
% Filename.
% The name of model file begins with an alphabetic character,
......@@ -21,24 +18,7 @@ function dynare(fname, varargin)
%
% part of DYNARE, copyright Dynare Team (2001-2008)
% Gnu Public License.
MATLAB = ver('matlab');
% FIXME:
% It's not satisfactory to convert string versions into numbers, and to
% compare these numbers:
% - conversion will fail if version = 1.2.3
% - it will give 7.10 < 7.9
VERSION = str2num(MATLAB.Version);
dynareroot = strrep(which('dynare.m'),'dynare.m','');
if (VERSION <= 7.4)
addpath([dynareroot '../mex/2007a/'])
else
addpath([dynareroot '../mex/2007b/'])
end
dynareroot = dynare_config();
if ~isstr(fname)
error ('The argument in DYNARE must be a text string.') ;
end
......
function dynareroot = dynare_config(path_to_dynare)
% This function tests the existence of valid mex files (for qz
% decomposition, solution to sylvester equation and kronecker
% products...) and, if needed, add paths to the matlab versions
% of these routines.
%
%
% INPUTS
% none
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
%
% part of DYNARE, copyright Dynare Team (2001-2008)
% Gnu Public License.
MATLAB = ver('matlab');
%% FIXME:
%% It's not satisfactory to convert string versions into numbers, and to
%% compare these numbers:
%% - conversion will fail if version = 1.2.3
%% - it will give 7.10 < 7.9
VERSION = str2num(MATLAB.Version);
if nargin
addpath(path_to_dynare);
end
dynareroot = strrep(which('dynare.m'),'dynare.m','');
if (VERSION <= 7.4)
addpath([dynareroot '../mex/2007a/'])
else
addpath([dynareroot '../mex/2007b/'])
end
%% Set mex routine names
mex_status = cell(1,3);
mex_status(1,1) = {'mjdgges'};
mex_status(1,2) = {'qz'};
mex_status(1,3) = {'Generalized QZ decomposition'};
mex_status(2,1) = {'gensylv'};
mex_status(2,2) = {'gensylv'};
mex_status(2,3) = {'Sylvester equation solution'};
mex_status(3,1) = {'A_times_B_kronecker_C'};
mex_status(3,2) = {'kronecker'};
mex_status(3,3) = {'Kronecker products'};
mex_status(4,1) = {'sparse_hessian_times_B_kronecker_C'};
mex_status(4,2) = {'kronecker'};
mex_status(4,3) = {'Sparse kronecker products'};
number_of_mex_files = size(mex_status,1);
%% Remove some directories from matlab's path. This is necessary if the user has
%% added dynare_v4/matlab with the subfolders. Matlab has to ignore these
%% subfolders if valid mex files exist.
matlab_path = path;
for i=1:number_of_mex_files
test = strfind(matlab_path,['/matlab/' mex_status{i,2}]);
action = length(test);
if action==1
rmpath([dynareroot mex_status{i,2}]);
matlab_path = path;
elseif action>1
disp('Dynare configuration problem!')
end
end
%% Test if valid mex files are available, if a mex file is not available
%% a matlab version of the routine is included in the path.
disp(' ')
disp('Configuring Dynare ...')
for i=1:number_of_mex_files
test = (exist(mex_status{i,1}) == 3);
if ~test
addpath([dynareroot mex_status{i,2}]);
message = '[m] ';
else
message = '[mex] ';
end
disp([ message mex_status{i,3} '.' ])
disp([])
end
disp(' ')
\ No newline at end of file
% solves a*x+b*x*c=d
function x=sylvester3(a,b,c,d)
n = size(a,1);
m = size(c,1);
if n == 1
x=d./(a*ones(1,m)+b*c);
return
end
if m == 1
x = (a+c*b)\d;
return;
end
x=zeros(n,m);
[u,t]=schur(c);
[aa,bb,qq,zz]=qz(full(a),full(b),'real'); % available in Matlab version 6.0
d=qq*d*u;
i = 1;
while i < m
if t(i+1,i) == 0
if i == 1
c = zeros(n,1);
else
c = bb*(x(:,1:i-1)*t(1:i-1,i));
end
x(:,i)=(aa+bb*t(i,i))\(d(:,i)-c);
i = i+1;
else
if i == n
c = zeros(n,1);
c1 = zeros(n,1);
else
c = bb*(x(:,1:i-1)*t(1:i-1,i));
c1 = bb*(x(:,1:i-1)*t(1:i-1,i+1));
end
z = [aa+bb*t(i,i) bb*t(i+1,i); bb*t(i,i+1) aa+bb*t(i+1,i+1)]...
\[d(:,i)-c;d(:,i+1)-c1];
x(:,i) = z(1:n);
x(:,i+1) = z(n+1:end);
i = i + 2;
end
end
if i == m
c = bb*(x(:,1:m-1)*t(1:m-1,m));
x(:,m)=(aa+bb*t(m,m))\(d(:,m)-c);
end
x=zz*x*u';
% 01/25/03 MJ corrected bug for i==m (sign of c in x determination)
% 01/31/03 MJ added 'real' to qz call
% solves iteratively ax+bxc=d
function x=sylvester3a(x0,a,b,c,d)
a_1 = inv(a);
b = a_1*b;
d = a_1*d;
e = 1;
iter = 1;
while e > 1e-8 & iter < 500
x = d-b*x0*c;
e = max(max(abs(x-x0)));
x0 = x;
iter = iter + 1;
end
if iter == 500
warning('sylvester3a : Only accuracy of %10.8f is achieved after 500 iterations')
end
\ No newline at end of file
......@@ -31,7 +31,7 @@ end
[mB,nB] = size(B);
if nargin == 3
[mC,nC] = size(C);
if mB*mc ~= nA
if mB*mC ~= nA
error('Input dimension error!')
end
D = zeros(mA,nB*nC);
......@@ -40,7 +40,7 @@ else
if mB*mB ~= nA
error('Input dimension error!')
end
D = D = zeros(mA,nB*nB);
D = zeros(mA,nB*nB);
loop = (mB*nB*mB*nB > 1e7);
end
% Computational part.
......
......@@ -18,5 +18,11 @@ function D = sparse_hessian_times_B_kronecker_C(A,B,C)
%
% part of DYNARE, copyright Dynare Team (1996-2008)
% Gnu Public License.
D = A_times_B_kronecker_C(A,B,C);
\ No newline at end of file
switch nargin
case 3
D = A_times_B_kronecker_C(A,B,C);
case 2
D = A_times_B_kronecker_C(A,B,B);
otherwise
error('Two or Three input arguments required!')
end
\ No newline at end of file
......@@ -34,7 +34,7 @@ end
% Check the first two inputs.
[me,ne] = size(e);
[md,nd] = size(d);
if ( ~isdouble(e) | ~isdouble(d) | iscomplex(e) | iscomplex(d) | me~=ne | md~=nd | me~=nd)
if ( ~isreal(e) | ~isreal(d) | me~=ne | md~=nd | me~=nd)
% info should be negative in this case, see dgges.f.
error('MYDGGES requires two square real matrices of the same dimension.')
end
......@@ -52,7 +52,7 @@ info = 0;
% Computational part.
try
[ss,tt,qq,w] = qz(e,d);
[ss,tt,qq,w] = qzdiv(qz_criterium,tt,ss,qq,w);
[tt,ss,qq,w] = qzdiv(qz_criterium,tt,ss,qq,w);
warning_old_state = warning;
warning off;
eigval = diag(ss)./diag(tt);
......
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