dynare_v4: adding scilab directory
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@809 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
michel
parent
d19eb2f7ce
commit
5edb7d1d00
|
|
@ -0,0 +1,2 @@
|
|||
function y=%s_2_s(a,b)
|
||||
y = real(a) > real(b)
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
copy a:*.bin c:\dynare\scilab
|
||||
copy a:*.sci c:\dynare\scilab
|
||||
copy a:*.exe c:\dynare\scilab
|
||||
copy a:names c:\dynare\scilab
|
||||
copy a:lib c:\dynare\scilab
|
||||
attrib -r c:\deamq\*
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
function [d1]=bksupk(ny,fid,jcf,icc1)
|
||||
d1=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('ykmax_','iter_','c','ncc');
|
||||
|
||||
icf = 1:jcf-1;
|
||||
ir = (iter_-1)*ny+1:ny*iter_;
|
||||
irf = icc1+(iter_-1)*ny;
|
||||
d1 = zeros(iter_*ny,1);
|
||||
|
||||
ofs = ((iter_-1)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v = c(:,jcf)
|
||||
d1(ir,1) = %v(:);
|
||||
ir = ir-ny;
|
||||
|
||||
i = 2;
|
||||
|
||||
while i<=ykmax_|i<=iter_ then
|
||||
irf1 = selif(irf,irf<iter_*ny);
|
||||
|
||||
ofs = ((iter_-i)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v1 = c(:,jcf)-c(:,1:size(irf1,1))*d1(irf1)
|
||||
d1(ir,1) = %v1(:);
|
||||
ir = ir-ny;
|
||||
irf = irf-ny;
|
||||
i = i+1;
|
||||
end
|
||||
|
||||
while i<=iter_ then
|
||||
|
||||
ofs = ((iter_-i)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v1 = c(:,jcf)-c(:,icf)*d1(irf)
|
||||
d1(ir,1) = %v1(:);
|
||||
ir = ir-ny;
|
||||
irf = irf-ny;
|
||||
i = i+1;
|
||||
end
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
function [d]=bksup1(ny,jcf)
|
||||
d=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global c
|
||||
|
||||
ir = (iter_-2)*ny+1:ny+(iter_-2)*ny;
|
||||
%irf = iyf+(iter_-1)*ny;
|
||||
icf = 1:size(iyf,2);
|
||||
|
||||
for i = 2:iter_
|
||||
c(ir,jcf) = c(ir,jcf)-c(ir,icf)*c(%irf,jcf);
|
||||
ir = ir-ny;
|
||||
%irf = %irf-ny;
|
||||
end
|
||||
|
||||
d = c(:,jcf);
|
||||
|
||||
|
|
@ -0,0 +1,53 @@
|
|||
function [d1]=bksupk(ny,fid,jcf,icc1)
|
||||
d1=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('ykmax_','iter_','c','ncc');
|
||||
|
||||
icf = 1:jcf-1;
|
||||
ir = (iter_-1)*ny+1:ny*iter_;
|
||||
%irf = icc1+(iter_-1)*ny;
|
||||
d1 = zeros(iter_*ny,1);
|
||||
|
||||
ofs = ((iter_-1)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v = c(:,jcf)
|
||||
d1(ir,1) = %v(:);
|
||||
ir = ir-ny;
|
||||
|
||||
i = 2;
|
||||
|
||||
while i<=ykmax_|i<=iter_ then
|
||||
irf1 = selif(%irf,%irf<iter_*ny);
|
||||
|
||||
ofs = ((iter_-i)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v1 = c(:,jcf)-c(:,1:size(irf1,1))*d1(irf1)
|
||||
d1(ir,1) = %v1(:);
|
||||
ir = ir-ny;
|
||||
%irf = %irf-ny;
|
||||
i = i+1;
|
||||
end
|
||||
|
||||
while i<=iter_ then
|
||||
|
||||
ofs = ((iter_-i)*ny+1-1)*ncc*8;
|
||||
junk = mseek(ofs,fid,'set');
|
||||
c = mtlb_fread(fid,[ncc,ny],'float64');
|
||||
c = c';
|
||||
|
||||
%v1 = c(:,jcf)-c(:,icf)*d1(%irf)
|
||||
d1(ir,1) = %v1(:);
|
||||
ir = ir-ny;
|
||||
%irf = %irf-ny;
|
||||
i = i+1;
|
||||
end
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
function [e_1,e_2,e_inf]=brm(v,nbr,drop)
|
||||
e_1=[];e_2=[];e_inf=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global('y_','iy_','iter_','ykmin_','ykmax_','it_','endo_nbr','ex_');
|
||||
|
||||
i = (iy_>0)';
|
||||
iyr0 = find(i(:))';
|
||||
|
||||
ex_old = ex_;
|
||||
|
||||
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
|
||||
ex_ = ex_(drop+1:mtlb(end),:);
|
||||
|
||||
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
|
||||
y = y_(:,drop+1:mtlb(end));
|
||||
y = y(:);
|
||||
|
||||
iter = iter_-drop-ykmin_-ykmax_;
|
||||
|
||||
z = zeros(endo_nbr,iter);
|
||||
|
||||
for it_ = ykmin_+1:iter+ykmin_
|
||||
z(:,it_) = ff_(y(iyr0));
|
||||
iyr0 = iyr0+endo_nbr;
|
||||
end
|
||||
|
||||
t1 = z(nbr,:)';
|
||||
|
||||
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
|
||||
t2 = v(drop+1:mtlb(end)-2);
|
||||
t = t1 ./ t2;
|
||||
|
||||
//! mtlb_mean(abs(t)) may be replaced by
|
||||
//! mean(abs(t)) if abs(t)is a vector
|
||||
//! mean(abs(t),1) if abs(t)is a matrix
|
||||
e_1 = log(mtlb_mean(abs(t)))/log(10);
|
||||
|
||||
//!! Unknown function var ,the original calling sequence is used
|
||||
e_2 = log(var(t))/log(10);
|
||||
|
||||
//! mtlb_max(abs(t)) may be replaced by
|
||||
//! max(abs(t)) if abs(t)is a vector
|
||||
//! max(abs(t),'r') if abs(t)is a matrix
|
||||
e_inf = log(mtlb_max(abs(t)))/log(10);
|
||||
|
||||
ex_ = ex_old;
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
function [x]=bseastr(s1,s2)
|
||||
x=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
m = size(s1,1);
|
||||
x = zeros(m,1);
|
||||
|
||||
//!! Unknown function deblank ,the original calling sequence is used
|
||||
s1 = convstr(deblank(s1),'u');
|
||||
|
||||
//!! Unknown function deblank ,the original calling sequence is used
|
||||
s2 = convstr(deblank(s2),'u');
|
||||
|
||||
for im = 1:m
|
||||
key = part(s1(im),:);
|
||||
h = size(s2,1);
|
||||
l = 1;
|
||||
while l<=h then
|
||||
mid = round((h+l)/2);
|
||||
temp = part(s2(mid),:);
|
||||
if ~(key==temp) then
|
||||
for i = 1:min(length(key),length(temp))
|
||||
if part(temp,i)>part(key,i) then
|
||||
h = mid-1;
|
||||
break
|
||||
;
|
||||
|
||||
else
|
||||
l = mid+1;
|
||||
break
|
||||
;
|
||||
|
||||
end
|
||||
end
|
||||
else
|
||||
x(im) = mid
|
||||
break
|
||||
;
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,188 @@
|
|||
function [result]=check()
|
||||
result=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
//
|
||||
|
||||
global ex_ it_ jacobia_
|
||||
|
||||
xlen = xkmax_+xkmin_+1;
|
||||
[klen,ny] = size(iy_);
|
||||
iyr0 = matrix(find(matrix(iy_',(ykmin_+ykmax_+1)*ny,1)),1,-1);
|
||||
it_ = ykmin_+1;
|
||||
tempex = ex_;
|
||||
|
||||
if valf_==0 then
|
||||
ys = ones(ykmin_+ykmax_+1,1).*.ys_;
|
||||
z = ys(iyr0);
|
||||
if exo_nbr > 0 then
|
||||
ex_ = ones(xlen,1)*exe_';
|
||||
end
|
||||
i = evstr(fname_+'_ff(z)');
|
||||
else
|
||||
t = y_(:,1:ykmin_+ykmax_+1);
|
||||
z = matrix(t,size(t,1)*size(t,2),1);
|
||||
z = z(iyr0);
|
||||
end
|
||||
|
||||
if ~(iy_(ykmin_+1,:)>0) then
|
||||
error('Model specification error: some variable doesn''t appear in current period.');
|
||||
|
||||
end
|
||||
|
||||
if (ykmax_==0)&(ykmin_==0) then
|
||||
error('Static model: no eigenvalue to compute.');
|
||||
end
|
||||
|
||||
it_ = ykmin_+1;
|
||||
jacob(fname_+'_ff',z);
|
||||
ex_ = tempex;
|
||||
clear('tempex');
|
||||
|
||||
k1 = iy_(find((((1:ykmin_+ykmax_+1)==(ykmin_+1))==%f)')',:);
|
||||
temp = matrix(k1',size(k1,1)*size(k1,2),1);
|
||||
k2 = temp(matrix(find(temp),1,-1));
|
||||
k3 = iy_(ykmin_+1,:)';
|
||||
a = jacobia_(:,k2);
|
||||
b = jacobia_(:,k3);
|
||||
a = b\a;
|
||||
clear('b','jacobia_');
|
||||
|
||||
if ykmin_>0 then
|
||||
k2 = cumsum(iy_(1:ykmin_,:),1);
|
||||
if ykmax_>0 then
|
||||
%v2 = iy_(ykmin_+2:klen,:)
|
||||
temp = %v2($:-1:1,:);
|
||||
if size(temp,1)>1 then
|
||||
temp = cumsum(temp,1);
|
||||
end
|
||||
k2 = [k2;temp($:-1:1,:)];
|
||||
end
|
||||
else
|
||||
%v1 = iy_(ykmin_+2:klen,:)
|
||||
temp = %v1($:-1:1,:);
|
||||
if size(temp,1)>1 then
|
||||
temp = cumsum(temp,1);
|
||||
end
|
||||
k2 = temp($:-1:1,:);
|
||||
end
|
||||
|
||||
%temp = k2';
|
||||
k3 = %temp(:);
|
||||
k4 = [find(k3)',k3(find(k3))];
|
||||
k3(k4(:,1)) = [1:size(k4,1)]';
|
||||
kad = zeros(0,0);
|
||||
|
||||
if ny<size(k3,1) then
|
||||
k4 = k3(1:size(k3,1)-ny,:)&k3(1+ny:size(k3,1),:);
|
||||
temp = [k3(1:size(k3,1)-ny,:),k4];
|
||||
kad = temp(find(temp(:,size(temp,2)))',:);
|
||||
if nnz(kad)==size(kad,1)*size(kad,2) then
|
||||
kad = kad(:,1);
|
||||
temp = [k3(1+ny:size(k3,1),:),k4];
|
||||
kae = temp(find(temp(:,size(temp,2)))',:);
|
||||
kae = kae(:,1);
|
||||
end
|
||||
end
|
||||
clear('temp','k4');
|
||||
|
||||
if ykmax_>0 then
|
||||
temp = [(1:ny)',k2([ykmin_,ykmin_+1],:)'];
|
||||
k3 = temp(matrix(find(sum(k2([ykmin_,ykmin_+1],:),1)),1,-1)',:);
|
||||
else
|
||||
temp = [(1:ny)',k2(ykmin_,:)'];
|
||||
k3 = temp(matrix(find(k2(ykmin_,:)'),1,-1),:);
|
||||
end
|
||||
|
||||
temp = matrix(k1',size(k1,1)*size(k1,2),1);
|
||||
temp = [temp,ones(klen-1,1).*.((1:ny)')];
|
||||
temp = [temp,((2:klen)').*.ones(ny,1)];
|
||||
k1 = temp(matrix(find(matrix(k2',size(k2,1)*size(k2,2),1)),1,-1),:);
|
||||
clear('temp');
|
||||
|
||||
nd = size(k1,1);
|
||||
d = zeros(nd,nd);
|
||||
e = d;
|
||||
|
||||
if ykmin_>0 then
|
||||
temp = [(1:nd)',k1(:,1)];
|
||||
temp1 = (k1(:,3)<=(ykmin_+1))&(k1(:,1)>0);
|
||||
k2 = temp(find(temp1)',:);
|
||||
e(1:size(k3,1),k2(:,1)) = -a(k3(:,1),k2(:,2));
|
||||
end
|
||||
|
||||
nfwrd = 0;
|
||||
if ykmax_>0 then
|
||||
temp = [(1:nd)',k1(:,1)];
|
||||
temp1 = (k1(:,3)>(ykmin_+1))&(k1(:,1)>0);
|
||||
k2 = temp(find(temp1)',:);
|
||||
d(1:size(k3,1),k2(:,1)) = a(k3(:,1),k2(:,2)-ny);
|
||||
nfwrd = sum(k1(:,3) > ykmin_+1);
|
||||
end
|
||||
|
||||
|
||||
if (ykmin_>0)&(ykmax_>0) then
|
||||
a = eye(size(k3,1),size(k3,1));
|
||||
temp = (1:nd)';
|
||||
temp1 = (1:size(k3,1))';
|
||||
k2 = [temp(find(k1(:,3)==(ykmin_+2))',:),temp1(find(k3(:,3))',:)];
|
||||
e(1:size(k3,1),k2(:,1)) = -a(:,k2(:,2));
|
||||
temp1 = [(1:size(k3,1))',k3(:,3)];
|
||||
k2 = [temp(find(k1(:,3)==(ykmin_+1))',:),temp1(find(k3(:,2))',:)];
|
||||
temp = k2(:,[1,2]);
|
||||
k2 = temp(find(k2(:,3)==0)',:);
|
||||
d(1:size(k3,1),k2(:,1)) = a(:,k2(:,2));
|
||||
clear('a');
|
||||
end
|
||||
|
||||
clear('temp','temp1');
|
||||
|
||||
if nnz(kad)==size(kad,1)*size(kad,2) then
|
||||
for j = 1:size(kad,1)
|
||||
d(size(k3,1)+j,kad(j)) = 1;
|
||||
e(size(k3,1)+j,kae(j)) = 1;
|
||||
end
|
||||
end
|
||||
|
||||
if grep(getversion(),'2.6') then
|
||||
[ss,tt] = gschur(e,d);
|
||||
else
|
||||
[ss,tt] = schur(e,d);
|
||||
end
|
||||
|
||||
ss = diag(ss);
|
||||
tt = diag(tt);
|
||||
|
||||
mod = ieee();
|
||||
ieee(2);
|
||||
|
||||
lambda = ss ./ tt;
|
||||
|
||||
ieee(mod);
|
||||
|
||||
%v = abs(lambda)
|
||||
[m_lambda,%i]=sort(%v);
|
||||
dyn_disp('Eigenvalues:');
|
||||
dyn_disp(' Modulus Real Imaginary');
|
||||
|
||||
mprintf('%12.6f %12.6f\n',[m_lambda lambda(%i)]);
|
||||
na = sum(abs(lambda)>1);
|
||||
mprintf('\nThere are %d eigenvalue(s) larger than 1 in modulus ',na);
|
||||
mprintf('for %d forward-looking variable(s)',nfwrd);
|
||||
mfprintf(fh_log,'%12.6f %12.6f\n',[m_lambda lambda(%i)]);
|
||||
mfprintf(fh_log,'\nThere are %d eigenvalue(s) larger than 1 in modulus ',na);
|
||||
mfprintf(fh_log,'for %d forward-looking variable(s)',nfwrd);
|
||||
|
||||
// 02/06/02 MJ port to Scilab
|
||||
|
||||
// 2/9/99 MJ: line 15, added test for absence of exogenous variable.
|
||||
// 8/27/2000 MJ: change JACOB call. Added ...,1 to cumsum()
|
||||
// 6/24/01 MJ: added count of abs(eigenvalues) > 1
|
||||
// 02/21/03 MJ: redone output to logfile
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
function []=d_corr(varargin)
|
||||
|
||||
if length(varargin) == 0
|
||||
for i=1:size(lgy_,1)
|
||||
varargin(i) = lgy_(i);
|
||||
end
|
||||
end
|
||||
|
||||
m = length(varargin);
|
||||
k = [];
|
||||
if m > 0 then
|
||||
dyn_disp(m)
|
||||
for i=1:m
|
||||
s = varargin(i)
|
||||
k = [k; grep_exact(lgy_,s)];
|
||||
end
|
||||
end
|
||||
|
||||
if k==[] then
|
||||
error('One of the variable specified does not exist');
|
||||
end
|
||||
|
||||
y = y_(k,:);
|
||||
y = y(:,100:$);
|
||||
n = size(y,2);
|
||||
y = y-mean(y,2)*ones(1,n);
|
||||
v = sum(y.*y,2);
|
||||
c = (y*y')./sqrt(v*v');
|
||||
|
||||
dyn_disp('Correlation matrix')
|
||||
mprintf('%8s ','')
|
||||
mfprintf(fh_log,'%8s ','')
|
||||
for i=1:m
|
||||
mprintf('%8s ',varargin(i))
|
||||
mfprintf(fh_log'%8s ',varargin(i))
|
||||
end
|
||||
mprintf('\n');
|
||||
mfprintf(fh_log,'\n');
|
||||
for i=1:m
|
||||
mprintf('%8s ',varargin(i))
|
||||
mfprintf(fh_log,'%8s ',varargin(i))
|
||||
for j=1:m
|
||||
mprintf('%8.4f ',c(i,j))
|
||||
mfprintf(fh_log,'%8.4f ',c(i,j))
|
||||
end
|
||||
mprintf('\n');
|
||||
mfprintf(fh_log,'\n');
|
||||
end
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
function []=dcompare(s1)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('dsmpl_','iter_');
|
||||
global('nvx','nvy','x','y','lag1');
|
||||
|
||||
ftest(s1,0);
|
||||
|
||||
i = (lag1(1):size(x,2)-lag1(2)+1)';
|
||||
|
||||
if size(dsmpl_,1)==1 then
|
||||
error('DSAMPLE not specified.');
|
||||
end
|
||||
|
||||
if dsmpl_(3)>0 then
|
||||
if dsmpl_(3)==2 then
|
||||
if dsmpl_(1)<0|dsmpl_(2)>(size(x,2)-lag1(2)) then
|
||||
error('Wrong sample.');
|
||||
end
|
||||
i = (dsmpl_(1)+lag1(1):dsmpl_(2)+lag1(1))';
|
||||
elseif dsmpl_(3)==1 then
|
||||
if dsmpl_(1)>(size(x,2)-lag1(2)) then
|
||||
error('Wrong sample.');
|
||||
end
|
||||
i = (lag1(1):dsmpl_(1)+lag1(1))';
|
||||
end
|
||||
end
|
||||
|
||||
j = bseastr(nvx,nvy);
|
||||
|
||||
|
||||
//! mtlb(stop) can be replaced by stop() or stop whether stop is an m-file or not
|
||||
if mtlb(stop) then
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
|
||||
//! mtlb_mean(abs(x(j,i)-y(j,i))) may be replaced by
|
||||
//! mean(abs(x(j,i)-y(j,i))) if abs(x(j,i)-y(j,i))is a vector
|
||||
//! mean(abs(x(j,i)-y(j,i)),1) if abs(x(j,i)-y(j,i))is a matrix
|
||||
z = mean(mtlb_mean(abs(x(j,i)-y(j,i))));
|
||||
|
||||
dyn_disp('The mean absolute difference between set '+s1(1,:)+'and set '+s1(2,:));
|
||||
dyn_disp('is : '+string(z));
|
||||
return
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,18 @@
|
|||
---------------------------------------------------------------------------
|
||||
DIFFEXT Numerical approximation for hessian.
|
||||
The method is Richardson`s extrapolation.
|
||||
Sample call
|
||||
[D,err,relerr,n] = diffext('f',x,delta,toler)
|
||||
Inputs
|
||||
f name of the function
|
||||
x differentiation point
|
||||
options matrix of algorithm parameters
|
||||
delta error goal (1e-12)
|
||||
toler relative error goal (1e-12)
|
||||
Return
|
||||
J Jacobian
|
||||
|
||||
NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
|
||||
|
||||
Modified F. Collard, August 2001
|
||||
---------------------------------------------------------------------------
|
||||
|
|
@ -0,0 +1,127 @@
|
|||
function [J]=diffext(f,x,options,varargin)
|
||||
[nargout,nargin] = argn(0)
|
||||
//---------------------------------------------------------------------------
|
||||
//DIFFEXT Numerical approximation for hessian.
|
||||
// The method is Richardson`s extrapolation.
|
||||
// Sample call
|
||||
// [D,err,relerr,n] = diffext('f',x,delta,toler)
|
||||
// Inputs
|
||||
// f name of the function
|
||||
// x differentiation point
|
||||
// options matrix of algorithm parameters
|
||||
// delta error goal (1e-12) (suppressed MJ 02/27/02)
|
||||
// toler relative error goal (1e-12)
|
||||
// Return
|
||||
// J Jacobian
|
||||
//
|
||||
// NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
|
||||
//
|
||||
// Modified F. Collard, August 2001
|
||||
//---------------------------------------------------------------------------
|
||||
if nargin > 2 then
|
||||
if ~(options==[]) then
|
||||
// %delta = options(1);
|
||||
toler = options(2);
|
||||
maxit = options(3);
|
||||
else
|
||||
// %delta = 1e-12;
|
||||
toler = 1e-12;
|
||||
maxit = 20;
|
||||
end
|
||||
else
|
||||
%delta = 1e-12;
|
||||
toler = 1e-12;
|
||||
maxit = 20;
|
||||
end
|
||||
con = 1.4;
|
||||
con2 = con*con;
|
||||
big = 1e30;
|
||||
safe = 2;
|
||||
if nargin > 3
|
||||
ff = evstr(f+'(x,varargin)');
|
||||
else
|
||||
ff = evstr(f+'(x)');
|
||||
end
|
||||
nx = size(x,1);
|
||||
nf = size(ff,1);
|
||||
J = zeros(nf,nx);
|
||||
for xi = 1:nx
|
||||
err = big*ones(nf,1);;
|
||||
// relerr = big*ones(nf,1);
|
||||
h = max(abs(x(xi))/10, 10*gstep_)*100*gstep_;
|
||||
dx = zeros(nx,1);
|
||||
dx(xi,1) = h;
|
||||
if nargin > 3
|
||||
fs = evstr(f+'(x+dx,varargin)');
|
||||
fm = evstr(f+'(x-dx,varargin)');
|
||||
else
|
||||
fs = evstr(f+'(x+dx)');
|
||||
fm = evstr(f+'(x-dx)');
|
||||
end
|
||||
D1 = (fs-fm)/(2*h);
|
||||
mask = ones(nf,1);
|
||||
j = 2;
|
||||
while (1)
|
||||
D2 = zeros(nf,j);
|
||||
h = h/con;
|
||||
dx(xi,1) = h;
|
||||
if nargin > 3
|
||||
fs = evstr(f+'(x+dx,varargin)');
|
||||
fm = evstr(f+'(x-dx,varargin)');
|
||||
else
|
||||
fs = evstr(f+'(x+dx)');
|
||||
fm = evstr(f+'(x-dx)');
|
||||
end
|
||||
for fi = 1:nf
|
||||
if mask(fi)
|
||||
err2 = big;
|
||||
D2(fi,1) = (fs(fi)-fm(fi))/(2*h);
|
||||
fac = con2;
|
||||
for k = 2:j
|
||||
D2(fi,k) = D2(fi,k-1)+(D2(fi,k-1)-D1(fi,k-1))/(fac-1);
|
||||
fac = con2*fac;
|
||||
errt = max(abs(D2(fi,k)-D2(fi,k-1)),abs(D2(fi,k)-D1(fi,k-1)));
|
||||
if errt <= err2
|
||||
err2 = errt;
|
||||
deriv = D2(fi,k);
|
||||
end
|
||||
end
|
||||
err(fi) = abs(D2(fi,j)-D1(fi,j-1));
|
||||
// relerr(fi) = 2*err(fi)/(abs(D2(fi,j))+abs(D1(fi,j-1))+%eps);
|
||||
// if (err(fi) < toler & relerr(fi) < %delta)| err(fi) > safe*err2
|
||||
if err(fi) < toler | err(fi) > safe*err2
|
||||
J(fi,xi) = deriv;
|
||||
mask(fi) = 0;
|
||||
end
|
||||
end
|
||||
end
|
||||
if (mask == 0) then
|
||||
break
|
||||
end
|
||||
j = j+1;
|
||||
if j == maxit
|
||||
error('DIFFEXT didn''t converge. Try to increase gstep_ (default 0.01)')
|
||||
end
|
||||
D1 = D2;
|
||||
end
|
||||
[m_err,i] = max(err);
|
||||
if m_err > 1e-12
|
||||
error('DIFFEXT obtains an accuracy > 1e-12. Try to increase gstep_ (default 0.01)')
|
||||
end
|
||||
// [m_err,i] = max(relerr);
|
||||
// if m_err > 1e-12
|
||||
// dyn_disp(err)
|
||||
// dyn_disp(relerr)
|
||||
// dyn_disp(D2)
|
||||
// error('DIFFEXT obtains an relative accuracy > 1e-12. Try to increase gstep_ (default 0.01)')
|
||||
// end
|
||||
end
|
||||
|
||||
// 10/12/2001 MJ modified initial h
|
||||
// 02/25/2002 MJ put equation look inside
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,200 @@
|
|||
function []=disp_dr(dr,order,var_list)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
nx = size(dr('ghx'),2);
|
||||
nu = size(dr('ghu'),2);
|
||||
|
||||
k1 = dr('order_var');
|
||||
kstate = dr('kstate');
|
||||
k = find(kstate(:,2) <= ykmin_+1);
|
||||
klag = kstate(k,[1 2]);
|
||||
|
||||
|
||||
nvar = size(var_list,1);
|
||||
if nvar == 0 then
|
||||
nvar = endo_nbr;
|
||||
ivar = [1:nvar];
|
||||
else
|
||||
ivar=zeros(nvar,1);
|
||||
for i=1:nvar
|
||||
i_tmp = grep_exact(lgy_(k1,:),var_list(i,:));
|
||||
if isempty(i_tmp) then
|
||||
dyn_disp(var_list(i,:));
|
||||
error (['One of the variable specified does not exist']) ;
|
||||
else
|
||||
ivar(i) = i_tmp;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
dyn_disp('POLICY AND TRANSITION FUNCTIONS');
|
||||
// variable names
|
||||
str = ' ';
|
||||
for i = 1:nvar
|
||||
str = str+msprintf('%16s',lgy_(k1(ivar(i))));
|
||||
end
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
//
|
||||
// constant
|
||||
//
|
||||
str = 'Constant ';
|
||||
flag = 0;
|
||||
for i = 1:nvar
|
||||
x = dr.ys(k1(ivar(i)));
|
||||
if order == 2 then
|
||||
x= x+dr.ghs2(ivar(i))/2;
|
||||
end
|
||||
if abs(x) > .000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
if order == 2 then
|
||||
str = '(Correction) ';
|
||||
flag = 0;
|
||||
for i = 1:nvar
|
||||
x = dr.ghs2(ivar(i))/2;
|
||||
if abs(x) > 1e-6 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x)
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
//
|
||||
// ghx
|
||||
//
|
||||
for k = 1:nx
|
||||
flag = 0;
|
||||
str1 = msprintf('%s(%d)',lgy_(k1(klag(k,1))),klag(k,2)-ykmin_-2);
|
||||
str = msprintf('%-20s',str1);
|
||||
for i = 1:nvar
|
||||
x = dr.ghx(ivar(i),k);
|
||||
if abs(x)>.000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
//
|
||||
// ghu
|
||||
//
|
||||
for k = 1:nu
|
||||
flag = 0;
|
||||
str = msprintf('%-20s',lgx_(k));
|
||||
for i = 1:nvar
|
||||
x = dr('ghu')(ivar(i),k);
|
||||
if abs(x)>.000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
|
||||
if order>1 then
|
||||
// ghxx
|
||||
for k = 1:nx
|
||||
for j = 1:k
|
||||
flag = 0;
|
||||
str = msprintf('%-20s',lgy_(k1(dr('nstatic')+k),:)+'(-1)'+lgy_(k1(dr('nstatic')+j),:)+'(-1)');
|
||||
for i = 1:endo_nbr
|
||||
if k==j then
|
||||
x = dr.ghxx(i,(k-1)*nx+j)/2;
|
||||
else
|
||||
x = dr.ghxx(i,(k-1)*nx+j);
|
||||
end
|
||||
if abs(x)>.000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
end
|
||||
//
|
||||
// ghuu
|
||||
//
|
||||
for k = 1:nu
|
||||
for j = 1:k
|
||||
flag = 0;
|
||||
str = msprintf('%-20s',lgx_(k,:)+' '+lgx_(j,:));
|
||||
for i = 1:endo_nbr
|
||||
if k==j then
|
||||
x = dr.ghuu(i,(k-1)*nu+j)/2;
|
||||
else
|
||||
x = dr.ghuu(i,(k-1)*nu+j);
|
||||
end
|
||||
if abs(x)>.000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
end
|
||||
//
|
||||
// ghxu
|
||||
//
|
||||
for k = 1:nx
|
||||
for j = 1:nu
|
||||
flag = 0;
|
||||
str = msprintf('%-20s',lgy_(k1(dr('nstatic')+k),:)+'(-1) '+lgx_(j,:));
|
||||
for i = 1:endo_nbr
|
||||
x = dr.ghxu(i,(k-1)*nu+j);
|
||||
if abs(x)>.000001 then
|
||||
flag = 1;
|
||||
str = str+msprintf('%16.6f',x);
|
||||
else
|
||||
str = str+' 0';
|
||||
end
|
||||
end
|
||||
if flag then
|
||||
mprintf('%s\n',str);
|
||||
mfprintf(fh_log,'%s\n',str);
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
// 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
|
||||
// 02/21/2003 MJ improved display format
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
function []=disp_moments(y,drop,var_list)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
nvar = size(var_list,1);
|
||||
if nvar == 0
|
||||
nvar = endo_nbr;
|
||||
ivar = [1:nvar]';
|
||||
else
|
||||
ivar=zeros(nvar,1);
|
||||
for i=1:nvar
|
||||
i_tmp = grep_exact(lgy_,var_list(i,:));
|
||||
if isempty(i_tmp)
|
||||
error (['One of the variable specified does not exist']) ;
|
||||
else
|
||||
ivar(i) = i_tmp;
|
||||
end
|
||||
end
|
||||
end
|
||||
y = y(ivar,drop+1:$)';
|
||||
|
||||
m = mean(y,1);
|
||||
y = y-ones(size(y,1),1).*.m;
|
||||
s2 = mean(y .* y,1);
|
||||
z = [m',s2',(mean(y.^3,1) ./ (s2.^1.5))',(mean(y.^4,1) ./ (s2 .* s2)-3)'];
|
||||
|
||||
dyn_disp('MOMENTS OF SIMULATED VARIABLES');
|
||||
dyn_disp('VARIABLE MEAN VARIANCE SKEWNESS KURTOSIS');
|
||||
for i = 1:nvar
|
||||
mprintf('%10s %10.6f %10.6f %10.6f %10.6f\n',lgy_(i,:),z(i,:))
|
||||
mfprintf(fh_log,'%10s %10.6f %10.6f %10.6f %10.6f\n',lgy_(i,:),z(i,:))
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
function []=disp_th_moments(dr,var_list)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global('endo_nbr','lgy_');
|
||||
|
||||
nvar = size(var_list,1);
|
||||
order = dr('order_var');
|
||||
if nvar==0 then
|
||||
nvar = endo_nbr;
|
||||
ivar = (1:nvar)';
|
||||
else
|
||||
ivar = zeros(nvar,1);
|
||||
for i = 1:nvar
|
||||
|
||||
//!! Unknown function strmatch ,the original calling sequence is used
|
||||
i_tmp = grep_exact(lgy_,var_list(i,:));
|
||||
if i_tmp==[] then
|
||||
error('One of the variable specified does not exist');
|
||||
else
|
||||
ivar(i,1) = i_tmp(:);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
m = dr('ys')(ivar);
|
||||
|
||||
|
||||
//!! Unknown function th_autocovariances ,the original calling sequence is used
|
||||
Gamma_y = th_autocovariances(dr,ivar);
|
||||
|
||||
s2 = diag(Gamma_y);
|
||||
z = [m,sqrt(s2),s2];
|
||||
|
||||
dyn_disp('THEORETICAL MOMENTS');
|
||||
dyn_disp('VARIABLE MEAN STD. DEV. VARIANCE');
|
||||
for i = 1:size(ivar,1)
|
||||
dyn_disp(sprintf('%16s %16.6f %16.6f %16.6f\n',lgy_(ivar(i),:),z(i,:)));
|
||||
end
|
||||
|
||||
// 10/09/02 MJ
|
||||
// 10/18/02 MJ added th_autocovariances() and provided for lags on several
|
||||
// periods
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,292 @@
|
|||
function [dr]=dr1(iorder,dr,correct)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
xlen = xkmax_+xkmin_+1;
|
||||
klen = ykmin_+ykmax_+1;
|
||||
iyv = iy_';
|
||||
iyv = iyv(:);
|
||||
iyr0 = find(iyv)';
|
||||
it_ = ykmin_+1;
|
||||
|
||||
|
||||
if exo_nbr==0 then
|
||||
exe_ = [];
|
||||
end
|
||||
|
||||
if ~(iy_(ykmin_+1,:)>0) then
|
||||
error('Error in model specification: some variables don""t appear as current');
|
||||
end
|
||||
|
||||
if ykmax_==0 then
|
||||
error('No forward variable: no point in using dr1');
|
||||
end
|
||||
|
||||
if xlen>1 then
|
||||
error('SS: stochastic exogenous variables must appear only at the'+' current period. Use additional endogenous variables');
|
||||
end
|
||||
|
||||
tempex = ex_;
|
||||
|
||||
it_ = ykmin_+1;
|
||||
|
||||
z = dr('ys')*ones(1,klen);
|
||||
z = z(iyr0);
|
||||
|
||||
jacobia_ = real(jacob_a('ff1_',[z;exe_]));
|
||||
|
||||
ex_ = tempex;
|
||||
tempex = [];
|
||||
|
||||
nz = size(z,1);
|
||||
pred_var = find(or(iy_(1:ykmin_,:),1));
|
||||
|
||||
fwrd_var = find(or(iy_(ykmin_+2:ykmin_+ykmax_+1,:),1));
|
||||
both_var = intersect(pred_var,fwrd_var);
|
||||
if size(both_var,'*')
|
||||
pred_var = setdiff(pred_var,both_var);
|
||||
fwrd_var = setdiff(fwrd_var,both_var);
|
||||
end
|
||||
nboth = length(both_var);
|
||||
npred = length(pred_var);
|
||||
nfwrd = length(fwrd_var);
|
||||
stat_var = setdiff(1:endo_nbr,union(union(pred_var,both_var),fwrd_var));
|
||||
nstatic = length(stat_var);
|
||||
order_var = [stat_var,pred_var,both_var,fwrd_var];
|
||||
k1 = iy_(find((1:klen)~=(ykmin_+1)),:);
|
||||
b = jacobia_(:,iy_(ykmin_+1,order_var));
|
||||
a = b\jacobia_(:,nonzeros(k1'));
|
||||
if exo_nbr then
|
||||
fu = b\jacobia_(:,nz+1:size(jacobia_,2));
|
||||
end
|
||||
|
||||
// building kmask for z state vector in t+1
|
||||
if ykmin_>0 then
|
||||
if ykmax_>0 then
|
||||
%v2 = iy_(ykmin_+2:size(iy_,1),order_var)
|
||||
kmask = cumsum(%v2($:-1:1,:),1);
|
||||
%v2 = cumsum(iy_(1:ykmin_,order_var),1)
|
||||
kmask = [kmask;%v2($:-1:1,:)];
|
||||
end
|
||||
else
|
||||
%v1 = iy_(ykmin_+1:klen,order_var)
|
||||
kmask = cumsum(%v1($:-1:1,:),1);
|
||||
end
|
||||
|
||||
kmask = kmask';
|
||||
kmask = kmask(:);
|
||||
i_kmask = find(kmask)';
|
||||
// index of nonzero entries in kmask
|
||||
nd = size(i_kmask,1);
|
||||
// size of the state vector
|
||||
%v = 1:nd
|
||||
kmask(i_kmask,1) = %v(:);
|
||||
|
||||
// auxiliary equations
|
||||
// elements that are both in z(t+1) and z(t)
|
||||
k1 = find(kmask(1:length(kmask)-endo_nbr)&kmask(endo_nbr+1:length(kmask)))';
|
||||
kad = [];
|
||||
kae = [];
|
||||
if ~(k1==[]) then
|
||||
kad = kmask(k1+endo_nbr);
|
||||
kae = kmask(k1);
|
||||
end
|
||||
|
||||
// composition of state vector
|
||||
// col 1: variable; col 2: lead/lag in z(t+1);
|
||||
// col 3: A cols for t+1 (D); col 4: A cols for t (E)
|
||||
|
||||
kstate = [ones(klen-1,1).*.[1:endo_nbr]' ((klen:-1:2)').*.ones(endo_nbr,1) zeros((klen-1)*endo_nbr,2)];
|
||||
%v = iy_(:,order_var)
|
||||
kiy = %v($:-1:1,:)';
|
||||
kiy = kiy(:);
|
||||
kstate(1:ykmax_*endo_nbr,3) = kiy(1:ykmax_*endo_nbr)-endo_nbr;
|
||||
kstate(ykmax_*endo_nbr+1:size(kstate,1),4) = kiy((ykmax_+1)*endo_nbr+1:length(kiy));
|
||||
// put in E only the current variables that are not already in D
|
||||
kstate = kstate(i_kmask,:);
|
||||
dr('kstate') = kstate;
|
||||
|
||||
sdyn = endo_nbr-nstatic;
|
||||
|
||||
// buildind D and E
|
||||
d = zeros(nd,nd);
|
||||
e = d;
|
||||
|
||||
k = find((kstate(:,2)>=(ykmin_+2)) & kstate(:,3))';
|
||||
d(1:sdyn,k) = a(nstatic+1:size(a,1),kstate(k,3));
|
||||
k1 = find(kstate(:,2)==(ykmin_+2))';
|
||||
|
||||
a1 = eye(sdyn,sdyn);
|
||||
e(1:sdyn,1:size(k1,1)) = -a1(:,kstate(k1,1)-nstatic);
|
||||
k = find((kstate(:,2)<=(ykmin_+1))&kstate(:,4))';
|
||||
e(1:sdyn,k) = -a(nstatic+1:size(a,1),kstate(k,4));
|
||||
k2 = find(kstate(:,2)==(ykmin_+1))';
|
||||
[k3,k4] = setdiff(kstate(k2,1),kstate(k1,1));
|
||||
|
||||
if size(k3,'*') > 0 then
|
||||
d(1:sdyn,k2(k4)) = a1(:,k3-nstatic);
|
||||
end
|
||||
|
||||
if ~(kad==[]) then
|
||||
for j = 1:size(kad,1)
|
||||
d(sdyn+j,kad(j)) = 1;
|
||||
e(sdyn+j,kae(j)) = 1;
|
||||
end
|
||||
end
|
||||
|
||||
if grep(getversion(),'2.6') then
|
||||
[ss,tt,w,sdim] = gschur(e,d,'d');
|
||||
else
|
||||
[ss,tt,w,sdim] = schur(e,d,'d');
|
||||
end
|
||||
|
||||
nba = nd-sdim;
|
||||
nyf = sum(kstate(:,2)>(ykmin_+1));
|
||||
if nba~=nyf then
|
||||
dyn_disp('LSS warning: Blanchard-Kahn conditions are not satisfied. Run CHECK to learn more!');
|
||||
halt();
|
||||
|
||||
end
|
||||
|
||||
np = nd-nyf;
|
||||
n2 = np+1;
|
||||
n3 = nyf;
|
||||
n4 = n3+1;
|
||||
// derivatives with respect to dynamic state variables
|
||||
// forward variables
|
||||
gx = -w(1:n3,n2:nd)'\w(n4:nd,n2:nd)';
|
||||
// predetermined variables
|
||||
hx = w(1:n3,1:np)'*gx+w(n4:nd,1:np)';
|
||||
hx = tt(1:np,1:np)*hx\(ss(1:np,1:np)*hx);
|
||||
|
||||
|
||||
k1 = find(kstate(n4:nd,2) == ykmin_+1);
|
||||
k2 = find(kstate(1:n3,2) == ykmin_+2);
|
||||
dr('ghx') = [hx(k1,:);gx(k2(nboth+1:$),:)]
|
||||
|
||||
// lead variables actually present in the model
|
||||
j3 = nonzeros(kstate(:,3));
|
||||
j4 = find(kstate(:,3));
|
||||
// derivatives with respect to exogenous variables
|
||||
if exo_nbr then
|
||||
a1 = eye(endo_nbr,endo_nbr);
|
||||
aa1 = [];
|
||||
if nstatic>0 then
|
||||
aa1 = a1(:,1:nstatic);
|
||||
end
|
||||
dr('ghu') = -[aa1,a(:,j3)*gx(j4,1:npred+nboth)+a1(:,nstatic+1:nstatic+npred+nboth),a1(:,nstatic+npred+nboth+1:$)]\fu
|
||||
end
|
||||
|
||||
// static variables
|
||||
if nstatic>0 then
|
||||
temp = -a(1:nstatic,j3)*gx(j4,:)*hx;
|
||||
j5 = find(kstate(n4:nd,4));
|
||||
temp(:,j5) = temp(:,j5)-a(1:nstatic,nonzeros(kstate(:,4)));
|
||||
dr.ghx = [temp; dr.ghx];
|
||||
temp = [];
|
||||
end
|
||||
|
||||
dr('order_var') = order_var
|
||||
dr('nstatic') = nstatic
|
||||
dr('npred') = npred+nboth
|
||||
|
||||
if iorder==1 then
|
||||
return
|
||||
end
|
||||
|
||||
// Second order
|
||||
tempex = ex_;
|
||||
|
||||
hessian = real(jacob2('ff1_',[z;exe_]));
|
||||
|
||||
ex_ = tempex;
|
||||
clear('tempex');
|
||||
|
||||
|
||||
np1 = sum(sum(iy_(1:ykmin_,:)>0,2));
|
||||
zx = [eye(npred+nboth,npred+nboth);dr('ghx');gx*hx;zeros(exo_nbr,np)];
|
||||
k1 = nonzeros(iy_(:,order_var)');
|
||||
if exo_nbr then
|
||||
k1 = [k1;(length(k1)+1:length(k1)+exo_nbr)'];
|
||||
end
|
||||
kk = matrix(1:length(k1)^2,length(k1),length(k1));
|
||||
// reorder hessian
|
||||
kk = kk(k1,k1);
|
||||
hessian = hessian(:,kk);
|
||||
rhs = -hessian*(zx.*.zx);
|
||||
rhs = rhs(:);
|
||||
|
||||
//lhs
|
||||
hxt = hx';
|
||||
gx = dr('ghx');
|
||||
gx = gx(size(gx,1)-nyf+1:size(gx,1),:);
|
||||
Inp = speye(np*np,np*np);
|
||||
|
||||
f1 = sparse(jacobia_(:,nonzeros(iy_(ykmin_+2,:))));
|
||||
// f_y'g_xx hx kron hx + f_y g_xx + f_y'g_x h_xx + f_x' h_xx
|
||||
//
|
||||
j1 = sparse(jacobia_(:,iy_(ykmin_+1,order_var)));
|
||||
lhs = Inp.*.j1;
|
||||
k = matrix(1:endo_nbr*np*np,endo_nbr,np*np);
|
||||
nr = size(k,1);
|
||||
lhs(:,k(nr-nyf+1:nr,:)) = lhs(:,k(nr-nyf+1:nr,:))+hxt.*.hxt.*.f1;
|
||||
lhs(:,k(nstatic+1:nstatic+np,:)) = lhs(:,k(nstatic+1:nstatic+np,:))+Inp.*.(f1*gx);
|
||||
dr('ghxx') = lhs\rhs
|
||||
dr('ghxx') = matrix(dr('ghxx'),endo_nbr,np*np)
|
||||
|
||||
//ghxu
|
||||
//rhs
|
||||
hu = dr('ghu')(nstatic+1:nstatic+np,:);
|
||||
zu = [zeros(np1,exo_nbr);dr('ghu');gx*hu;eye(exo_nbr,exo_nbr)];
|
||||
nr = size(dr('ghxx'),1);
|
||||
rhs = -hessian*(zx.*.zu)-f1*dr('ghxx')(nr-nyf+1:nr,:)*(hx.*.hu);
|
||||
|
||||
//lhs
|
||||
lhs = j1;
|
||||
lhs(:,nstatic+1:nstatic+np) = lhs(:,nstatic+1:nstatic+np)+f1*gx;
|
||||
dr('ghxu') = lhs\rhs
|
||||
|
||||
//ghuu
|
||||
//rhs
|
||||
rhs = -hessian*(zu.*.zu)-f1*dr('ghxx')(nr-nyf+1:nr,:)*(hu.*.hu);
|
||||
|
||||
//lhs
|
||||
dr('ghuu') = lhs\rhs
|
||||
|
||||
|
||||
// dr.ghs2
|
||||
// derivatives of F with respect to forward variables
|
||||
k1 = matrix(1:size(hessian,2),size(jacobia_,2),size(jacobia_,2));
|
||||
k1 = k1(npred+nboth+endo_nbr+1:npred+endo_nbr+2*nboth+nfwrd,npred+nboth+endo_nbr+1:npred+endo_nbr+2*nboth+nfwrd);
|
||||
|
||||
nr = size(dr('ghuu'),1);
|
||||
guu = dr('ghuu')(nr-nyf+1:nr,:);
|
||||
nr = size(dr('ghu'),1);
|
||||
gu = dr('ghu')(nr-nyf+1:nr,:);
|
||||
rhs = -(f1*guu+hessian(:,k1)*(gu.*.gu))*Sigma_e_(:);
|
||||
lhs = j1;
|
||||
lhs(:,nstatic+1:nstatic+np) = lhs(:,nstatic+1:nstatic+np)+f1*gx;
|
||||
nc = size(lhs,2);
|
||||
lhs(:,nc-nyf+1:nc) = lhs(:,nc-nyf+1:nc)+f1;
|
||||
dr('ghs2') = lhs\rhs;
|
||||
|
||||
// 01/08/2001 MJ put return on iorder == 1 after defining dr.kstate and dr.kdyn
|
||||
// 01/17/2001 MJ added dr.delta_s: correction factor for order = 2
|
||||
// 01/21/2001 FC correction of delta_s for more than 1 shock
|
||||
// 01/23/2001 MJ suppression of redundant sum() in delta_s formula
|
||||
// 02/22/2001 MJ stderr_ replaced by Sigma_e_
|
||||
// 08/24/2001 MJ changed the order of the variables, separates static
|
||||
// variables and handles only one instance of variables both
|
||||
// in lead and lag
|
||||
// 08/24/2001 MJ added sigma to state variables as in Schmitt-Grohe and
|
||||
// Uribe (2001)
|
||||
// 10/30/2002 MJ corrected bugs with lags on several periods
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
DSAMPLE : DSAMPLE(d1,d2)
|
||||
This optional command permits to reduce the number of
|
||||
periods considered in following output commands.
|
||||
If only one argument is
|
||||
provided, output is from period 1 to the period
|
||||
specified in the DSAMPLE command. If two arguments are
|
||||
present output is done for the interval between the
|
||||
two periods.
|
||||
DSAMPLE without arguments reset the sample to the one
|
||||
specified by PERIODS
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
function []=dsample(s1,s2)
|
||||
[nargout,nargin] = argn(0)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// DSAMPLE : DSAMPLE(d1,d2)
|
||||
// This optional command permits to reduce the number of
|
||||
// periods considered in following output commands.
|
||||
// If only one argument is
|
||||
// provided, output is from period -maxlag+1 to the period
|
||||
// specified in the DSAMPLE command. If two arguments are
|
||||
// present output is done for the interval between the
|
||||
// two periods.
|
||||
// DSAMPLE without arguments reset the sample to the one
|
||||
// specified by PERIODS
|
||||
|
||||
global dsmpl_
|
||||
|
||||
dsmpl_ = zeros(2,1);
|
||||
|
||||
if nargin==0 then
|
||||
dsmpl_(1) = 1
|
||||
dsmpl_(2,1) = iter_(:)+ykmin_+ykmax_;
|
||||
elseif nargin==1 then
|
||||
dsmpl_(1) = 1
|
||||
dsmpl_(2,1) = s1(:)+ykmin_;
|
||||
else
|
||||
dsmpl_(1,1) = s1(:)+ykmin_;
|
||||
dsmpl_(2,1) = s2(:)+ykmin_;
|
||||
end
|
||||
|
||||
if dsmpl_(1) > iter_+ykmin_|dsmpl_(2) > iter_+ykmin_+ykmax_ then
|
||||
t = 'DYNARE dsample error: one of the arguments is larger than the one'+' specified in PERIODS';
|
||||
error(t);
|
||||
end
|
||||
|
||||
// 02/23/01 MJ added error checking
|
||||
// 02/19/01 MJ added ykmin_
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
function dtranslatepaths(Paths,res_path)
|
||||
// Copyright INRIA
|
||||
if exists('m2scilib')==0 then load('SCI/macros/m2sci/lib'),end
|
||||
//logfile=%io(2)
|
||||
|
||||
|
||||
res_path=stripblanks(res_path)
|
||||
if part(res_path,length(res_path))<>'/' then
|
||||
res_path=res_path+'/'
|
||||
end
|
||||
Paths=stripblanks(Paths)
|
||||
logfile=file('open',res_path+'log','unknown')
|
||||
whsfil_unit=file('open',res_path+'whatis','unknown')
|
||||
|
||||
for k=1:size(Paths,'*')
|
||||
if part(Paths(k),length(Paths(k)))<>'/' then
|
||||
if MSDOS then
|
||||
Paths(k)=Paths(k)+'\',
|
||||
else
|
||||
Paths(k)=Paths(k)+'/',
|
||||
end
|
||||
end
|
||||
end
|
||||
mfiles=[]
|
||||
for k=1:size(Paths,'*')
|
||||
path=Paths(k)
|
||||
if MSDOS then
|
||||
mfiles=[mfiles;unix_g('dir /b '+path+'*.m')]
|
||||
sep='\'
|
||||
else
|
||||
mfiles=[mfiles;unix_g('ls '+path+'*.m')]
|
||||
sep='/'
|
||||
end
|
||||
end
|
||||
for k1=1:size(mfiles,1)
|
||||
if MSDOS then
|
||||
fnam=part(mfiles(k1),1:length(mfiles(k1))-2)
|
||||
mpath=path+mfiles(k1)
|
||||
else
|
||||
kk=strindex(mfiles(k1),sep)
|
||||
fnam=part(mfiles(k1),kk($)+1:length(mfiles(k1))-2)
|
||||
mpath=mfiles(k1)
|
||||
end
|
||||
scipath=res_path+fnam+'.sci'
|
||||
scepath=res_path+fnam+'.sce'
|
||||
if newest(mpath,scipath,scepath)==1 then
|
||||
mfile2sci(mpath,res_path,%f,%t)
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,58 @@
|
|||
function [z,zss]=dyn2vec(s1,s2)
|
||||
z=[];zss=[];
|
||||
[nargout,nargin] = argn(0)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
if dsmpl_==0 then
|
||||
%kk = 1:size(y_,2);
|
||||
else
|
||||
%kk = ykmin_+dsmpl_(1):ykmin_+dsmpl_(2);
|
||||
end
|
||||
|
||||
if nargin==0 then
|
||||
if nargout>1 then
|
||||
t = 'DYNARE dyn2vec error: the function doesn''t return values when'+' used without input argument';
|
||||
error(t);
|
||||
end
|
||||
for %i = 1:size(y_,1)
|
||||
tt = ['global '+string(lgy_(%i,:)) ; string(lgy_(%i,:))+' = y_(%i,%kk);'];
|
||||
deff('[]=assignin()',tt);
|
||||
assignin();
|
||||
clear('assignin')
|
||||
end
|
||||
return
|
||||
else
|
||||
j = grep_exact(lgy_,s1);
|
||||
if ~(j==[]) then
|
||||
%z = y_(j,%kk)';
|
||||
else
|
||||
j = grep_exact(lgx_,s1);
|
||||
if ~(j==[]) then
|
||||
if dsmpl_==0 then
|
||||
%z = ex_(:,j);
|
||||
else
|
||||
%z = ex_(xkmin_+dsmpl_(1):xkmin_+dsmpl_(2));
|
||||
end
|
||||
else
|
||||
t = 'DYNARE dyn2vec error: variable '+string(s1(%i,:))+' doesn''t'+' exist.';
|
||||
error(t);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
if nargout==1 then
|
||||
deff('[]=assignin()',['global '+string(s1) ; string(s1)+' = %z;']);
|
||||
assignin();
|
||||
clear('assignin')
|
||||
else
|
||||
zss = ys_(j);
|
||||
end
|
||||
|
||||
// 02/23/01 MJ redone, incorporating FC's improvements
|
||||
// 08/24/01 MJ replaced globalize by internal assignin
|
||||
// 08/24/01 MJ added 'exact' to strmatch (thanks to David Vavra)
|
||||
// 09/24/01 MJ translated to SciLab
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
function []=dyn_disp(s)
|
||||
disp(s)
|
||||
mputl(s,fh_log);
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
function []=dyn_mprintf(fmt,varargin)
|
||||
disp(fmt)
|
||||
disp(varargin)
|
||||
mprintf(fmt,varargin);
|
||||
mfprintf(fh_log,fmt,varargin);
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
DYNARE ( 'Filename' )
|
||||
This command runs the .MOD file specified in Filename.
|
||||
Filename could be enter with or without the .MOD extension.
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
function []=dynare(fname)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// DYNARE ( 'Filename' )
|
||||
// This command runs the .MOD file specified in Filename.
|
||||
// Filename could be enter with or without the .MOD extension.
|
||||
|
||||
if ~(type(fname)==10) then
|
||||
error('The argument in DYNARE must be a text string.');
|
||||
end
|
||||
k = strindex(fname,'.')
|
||||
if k then
|
||||
fname = part(fname,1:k-1);
|
||||
end
|
||||
|
||||
command = 'e:/dynare/scilab/dynare_s '+fname;
|
||||
|
||||
disp(command);
|
||||
|
||||
unix_w(command);
|
||||
|
||||
exec(fname+'.sci',-1) ;
|
||||
|
||||
endfunction
|
||||
// MJ 2/9/99: replace clear function by clear ff_
|
||||
// MJ 4/7/00: change the path of dynare_m
|
||||
// MJ 02/26/01: replaced local variable x by fname
|
||||
// MJ 09/14/01: Scilab translation
|
||||
Binary file not shown.
|
|
@ -0,0 +1,6 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
DYNASAVE : DYNASAVE ( [ 'filename' ] )
|
||||
This optional command saves the simulation results
|
||||
in a .BIN datafile. This command must follow SIMUL.
|
||||
Filename must be specified without the .BIN extension.
|
||||
|
|
@ -0,0 +1,36 @@
|
|||
function []=dynasave(s)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// DYNASAVE : DYNASAVE ( [ 'filename' ] )
|
||||
// This optional command saves the simulation results
|
||||
// in a .BIN datafile. This command must follow SIMUL.
|
||||
// Filename must be specified without the .BIN extension.
|
||||
|
||||
|
||||
if ~(matrix(find(abs(s)==46),1,-1)==[]) then
|
||||
dyn_disp('Warning : Error in DYNASAVE');
|
||||
dyn_disp('Filename '+s+' is incorrect.');
|
||||
dyn_disp('You must enter filename without extension,');
|
||||
dyn_disp('For example : '+s(1,1:matrix(find(abs(s)==46),1,-1)-1));
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
s = s+'.BIN';
|
||||
|
||||
[fid,%v] = mopen(s,'w',0)
|
||||
if %v<0 then fid = -1;end
|
||||
mtlb_fwrite(fid,size(y_),'int');
|
||||
mtlb_fwrite(fid,size(lgy_,2),'int');
|
||||
mtlb_fwrite(fid,ykmin_,'int');
|
||||
mtlb_fwrite(fid,ykmax_,'int');
|
||||
mtlb_fwrite(fid,xkmin_,'int');
|
||||
mtlb_fwrite(fid,xkmax_,'int');
|
||||
mtlb_fwrite(fid,lgy_,'int');
|
||||
mtlb_fwrite(fid,y_,'float64');
|
||||
mclose(fid);
|
||||
|
||||
return
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
DYNATYPE : DYNATYPE ( [ 'filename' ] )
|
||||
This optional command saves the simulation
|
||||
results in a text file. The name of each
|
||||
variable preceeds the corresponding results.
|
||||
This command must follow SIMUL.
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
function []=dynatype(s)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// DYNATYPE : DYNATYPE ( [ 'filename' ] )
|
||||
// This optional command saves the simulation
|
||||
// results in a text file. The name of each
|
||||
// variable preceeds the corresponding results.
|
||||
// This command must follow SIMUL.
|
||||
|
||||
[fid,%v] = mopen(s,'w',0)
|
||||
if %v<0 then fid = -1;end
|
||||
for i = 1:size(lgy_,1)
|
||||
mtlb_fprintf(fid,lgy_(i,:),'\n');
|
||||
mtlb_fprintf(fid,'\n');
|
||||
mtlb_fprintf(fid,'%15.8g\n',y_(i,:)');
|
||||
end
|
||||
mclose(fid);
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,11 @@
|
|||
function [y]=equiv(x)
|
||||
y=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
ys_ = x;
|
||||
|
||||
//!! Unknown function taylor ,the original calling sequence is used
|
||||
taylor();
|
||||
|
||||
//!! Unknown function newt ,the original calling sequence is used
|
||||
y = newt('ff20_',ys_);
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,11 @@
|
|||
function [y]=ff1_(x)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global ex_
|
||||
|
||||
n1 = length(x)-exo_nbr;
|
||||
ex_(it_+xkmin_-ykmin_,:) = x(n1+1:$)';
|
||||
y = evstr(fname_+'_ff(x(1:n1))');
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
function [y]=ff2_(x)
|
||||
y=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global('np_','endo_nbr','nf_','ex_','it_','ykmin_','xkmin_','exo_nbr','jp_','jf_','kf_','exe_');
|
||||
y = zeros(np_+endo_nbr+nf_,1);
|
||||
%v = x(1:np_)
|
||||
y(1:np_,1) = %v(:);
|
||||
ex_(it_-ykmin_+xkmin_,:) = x(np_+1:np_+exo_nbr)';
|
||||
|
||||
//!! Unknown function h1_ ,the original calling sequence is used
|
||||
%v = h1_(y(1:np_),ex_(it_-ykmin_+xkmin_,:)')
|
||||
y(np_+jp_,1) = %v(:);
|
||||
|
||||
//!! Unknown function g1_ ,the original calling sequence is used
|
||||
%v = g1_(y(1:np_),ex_(it_-ykmin_+xkmin_,:)')
|
||||
y(np_+jf_,1) = %v(:);
|
||||
|
||||
//!! Unknown function g1_ ,the original calling sequence is used
|
||||
%v = g1_(y(np_+jp_),exe_)
|
||||
y(kf_,1) = %v(:);
|
||||
y = ff_(y);
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
function [z]=ff2a_(y,dr)
|
||||
z=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
z = fff_(y)+dr('delta_s');
|
||||
|
|
@ -0,0 +1 @@
|
|||
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
function [z]=ff_(y)
|
||||
z=[];
|
||||
//
|
||||
z = zeros(2,1);
|
||||
|
||||
z(1) = y(2)+y(3)-aa*ex_(it_-1,1)*(y(1)^alph)-(1-delt)*y(1)
|
||||
%v = y(2)^(-gam)-((1+bet)^(-1))*(aa*alph*ex_(it_,1)*(y(3)^(alph-1))+1-delt)*(y(4)^(-gam))
|
||||
z(2,1) = %v(:);
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
function [x]=fff(y)
|
||||
x=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('ykmin_','ykmax_','iy_');
|
||||
|
||||
iyr = matrix(find(matrix(iy_'>0,(ykmin_+ykmax_+1)*size(iy_,2),1)),1,-1);
|
||||
y = ones(ykmin_+ykmax_+1,1).*.y;
|
||||
y = y(iyr,:);
|
||||
x = ff_(y);
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1 @@
|
|||
|
||||
|
|
@ -0,0 +1,10 @@
|
|||
function [z]=fff_(y)
|
||||
z=[];
|
||||
//
|
||||
z = zeros(2,1);
|
||||
global('ex_','it_','recur_');
|
||||
|
||||
global('aa','alph','bet','delt','gam');
|
||||
z(1) = y(1)+y(2)-aa*ex_(it_-1,1)*(y(2)^alph)-(1-delt)*y(2)
|
||||
%v = y(1)^(-gam)-((1+bet)^(-1))*(aa*alph*ex_(it_,1)*(y(2)^(alph-1))+1-delt)*(y(1)^(-gam))
|
||||
z(2,1) = %v(:);
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,24 @@
|
|||
function [a,b]=ffill(x,ixc,y)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
xc = size(x,1);
|
||||
|
||||
if y==[] then
|
||||
b = [ixc;0];
|
||||
a = [x,zeros(size(x,1),1)];
|
||||
else
|
||||
yc = size(y,1);
|
||||
b = [ixc;yc];
|
||||
|
||||
if xc>yc then
|
||||
a = [x,[y;zeros(xc-yc,size(y,2))]];
|
||||
elseif yc>xc then
|
||||
a = [[x;zeros(yc-xc,size(x,2))],y];
|
||||
else
|
||||
a = [x,y];
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
// 2001/09/1 MJ corrected for absent lags
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,73 @@
|
|||
function []=ftest(s1,s2)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('nvx','nvy','x','y','lag1');
|
||||
|
||||
if size(s1,1)~=2 then
|
||||
error('Spécifiez deux fichiers pour la comparaison.');
|
||||
end
|
||||
|
||||
for i = 1:2
|
||||
if ~(matrix(find(abs(s1(i,:))==46),1,-1)==[]) then
|
||||
error('Entrez les noms de fichiers sans extensions.');
|
||||
end
|
||||
end
|
||||
|
||||
s1 = s1+[' ';' '];
|
||||
file1 = part(s1(1),1:min(find(abs(part(s1(1),:))==32))-1)+'.BIN';
|
||||
file2 = part(s1(2),1:min(find(abs(part(s1(2),:))==32))-1)+'.BIN';
|
||||
|
||||
[fid,%v] = mopen(file1,'r',0)
|
||||
if %v<0 then fid = -1;end
|
||||
n1 = mtlb_fread(fid,1,'int');
|
||||
n2 = mtlb_fread(fid,1,'int');
|
||||
n3 = mtlb_fread(fid,1,'int');
|
||||
lag1 = mtlb_fread(fid,4,'int');
|
||||
nvx = mtlb_fread(fid,[n1,n3],'int');
|
||||
x = mtlb_fread(fid,[n1,n2],'float64');
|
||||
mclose(fid);
|
||||
nvx = ascii(nvx);
|
||||
|
||||
[fid,%v] = mopen(file2,'r',0)
|
||||
if %v<0 then fid = -1;end
|
||||
n1 = mtlb_fread(fid,1,'int');
|
||||
n2 = mtlb_fread(fid,1,'int');
|
||||
n3 = mtlb_fread(fid,1,'int');
|
||||
lag2 = mtlb_fread(fid,4,'int');
|
||||
nvy = mtlb_fread(fid,[n1,n3],'int');
|
||||
y = mtlb_fread(fid,[n1,n2],'float64');
|
||||
mclose(fid);
|
||||
nvy = ascii(nvy);
|
||||
|
||||
if size(x,1)~=size(y,1) then
|
||||
error('FTEST: The two files don''t have the same number of variables.');
|
||||
end
|
||||
|
||||
for i = 1:size(x,1)
|
||||
if ~(part(nvx(i),:)==part(nvy(i),:)) then
|
||||
error('FTEST: The two files don''t have the same variables.');
|
||||
end
|
||||
end
|
||||
|
||||
if nnz(lag1-lag2)>0 then
|
||||
error('FTEST: Leads and lags aren''t the same in both files.');
|
||||
end
|
||||
|
||||
j = zeros(size(s2,1),1);
|
||||
for i = 1:size(s2,1)
|
||||
|
||||
//!! Unknown function strmatch ,the original calling sequence is used
|
||||
k = strmatch(s2(i,:),nvx,'exact');
|
||||
if k==[] then
|
||||
t = 'FTEST: Variable '+s2(i)+'doesn''t exist';
|
||||
error(t);
|
||||
else
|
||||
j(i,1) = k(:);
|
||||
end
|
||||
end
|
||||
|
||||
y = y(j,:);
|
||||
x = x(j,:);
|
||||
|
||||
//06/18/01 MJ replaced beastr by strmatch
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
function [y]=fx_(x)
|
||||
y=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global('ys_','ykmin_','ykmax_','it_','iy_','xkmin_','ex_');
|
||||
|
||||
//!! Unknown function repmat ,the original calling sequence is used
|
||||
y = repmat(ys_,ykmin_+ykmax_+1,1);
|
||||
iyr0 = matrix(find(matrix(iy_',size(iy_,1)*size(iy_,2),1)),1,-1);
|
||||
y = y(iyr0);
|
||||
ex_(it_+xkmin_-ykmin_,:) = x';
|
||||
y = ff_(y);
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
GCOMPARE : GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
|
||||
This optional command plots the trajectories of a list of
|
||||
variables in two different simulations. One plot is drawn
|
||||
for each variable. The trajectories must have been previously
|
||||
saved by the instruction DYNASAVE. The simulation in file1
|
||||
is refered to as the base simulation.
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
function []=gcompare(s1,s2)
|
||||
// GCOMPARE : GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// This optional command plots the trajectories of a list of
|
||||
// variables in two different simulations. One plot is drawn
|
||||
// for each variable. The trajectories must have been previously
|
||||
// saved by the instruction DYNASAVE. The simulation in file1
|
||||
// is refered to as the base simulation.
|
||||
|
||||
|
||||
global('dsmpl_','iter_','ykmin_');
|
||||
global('nvx','nvy','x','y','lag1');
|
||||
|
||||
ftest(s1,s2);
|
||||
|
||||
ix = (1-lag1(1):size(x,2)-lag1(1))';
|
||||
i = (lag1(1):size(ix,1)-lag1(2)+1)';
|
||||
|
||||
if dsmpl_==0 then
|
||||
i = (ykmin_:size(y,2))';
|
||||
else
|
||||
i = (dsmpl_(1):dsmpl_(2))';
|
||||
end
|
||||
|
||||
for k = 1:size(x,1)
|
||||
xset('window',max(winsid()+1))
|
||||
max(winsid());
|
||||
mtlb_plot(ix(i),x(k,i),ix(i),y(k,i));
|
||||
xtitle(' ','Periods',' ');
|
||||
xtitle('Variable '+s2(k,:),' ',' ');
|
||||
l = min(i)+1;
|
||||
ll = max(i)-1;
|
||||
xstring(l,x(k,l),s1(1,:));
|
||||
xstring(ll,y(k,ll),s1(2,:));
|
||||
end
|
||||
|
||||
// 06/18/01 MJ corrected treatment of dsmpl_
|
||||
// 06/24/01 MJ removed color specification
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
function [row]=grep_exact(str1,str2)
|
||||
row = [];
|
||||
r = grep(str1,str2);
|
||||
for i=1:length(r)
|
||||
if str1(r(i)) == str2 then
|
||||
row = r(i);
|
||||
break
|
||||
end
|
||||
end
|
||||
|
|
@ -0,0 +1,18 @@
|
|||
---------------------------------------------------------------------------
|
||||
HESSEXT Numerical approximation for hessian.
|
||||
The method is Richardson`s extrapolation.
|
||||
Sample call
|
||||
[H] = hessext(f,x,options,varargin)
|
||||
Inputs
|
||||
f name of the function
|
||||
x differentiation point
|
||||
options matrix of algorithm parameters
|
||||
delta error goal (1e-12)
|
||||
toler relative error goal (1e-12)
|
||||
Return
|
||||
J Jacobian
|
||||
|
||||
NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
|
||||
|
||||
Modified F. Collard, August 2001
|
||||
---------------------------------------------------------------------------
|
||||
|
|
@ -0,0 +1,139 @@
|
|||
function [H]=hessext(f,x,options,varargin)
|
||||
[nargout,nargin] = argn(0)
|
||||
//---------------------------------------------------------------------------
|
||||
//HESSEXT Numerical approximation for hessian.
|
||||
// The method is Richardson`s extrapolation.
|
||||
// Sample call
|
||||
// [H] = hessext(f,x,options,varargin)
|
||||
// Inputs
|
||||
// f name of the function
|
||||
// x differentiation point
|
||||
// options matrix of algorithm parameters
|
||||
// %delta error goal (1e-12)
|
||||
// toler relative error goal (1e-12)
|
||||
// Return
|
||||
// J Jacobian
|
||||
//
|
||||
// NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
|
||||
//
|
||||
// Modified F. Collard, August 2001
|
||||
//---------------------------------------------------------------------------
|
||||
if nargin > 2 then
|
||||
if ~(options==[]) then
|
||||
%delta = options(1);
|
||||
toler = options(2);
|
||||
maxit = options(3)
|
||||
else
|
||||
%delta = 1e-9;
|
||||
toler = 1e-9;
|
||||
maxit = 20;
|
||||
end
|
||||
else
|
||||
%delta = 1e-9;
|
||||
toler = 1e-9;
|
||||
maxit = 20;
|
||||
end
|
||||
con = 1.4;
|
||||
con2 = con*con;
|
||||
big = 1e30;
|
||||
safe = 2;
|
||||
if nargin > 3
|
||||
ff = evstr(f+'(x,varargin)');
|
||||
else
|
||||
ff = evstr(f+'(x)');
|
||||
end
|
||||
nx = size(x,1);
|
||||
nf = size(ff,1);
|
||||
H = sparse([],[],[nf,nx*nx]);
|
||||
for xi = 1:nx
|
||||
for xj = xi:nx
|
||||
err = big*ones(nf,1);
|
||||
relerr = big*ones(nf,1);
|
||||
err2 = big;
|
||||
hx = max(abs(x(xi))/5,gstep_);
|
||||
hy = max(abs(x(xj))/5,gstep_);
|
||||
dx = zeros(nx,1);
|
||||
dy = zeros(nx,1);
|
||||
dx(xi) = hx;
|
||||
dy(xj) = hy;
|
||||
if nargin > 3
|
||||
fss = evstr(f+'(x+dx+dy,varargin)');
|
||||
fsm = evstr(f+'(x+dx-dy,varargin)');
|
||||
fms = evstr(f+'(x-dx+dy,varargin)');
|
||||
fmm = evstr(f+'(x-dx-dy,varargin)');
|
||||
else
|
||||
fss = evstr(f+'(x+dx+dy)');
|
||||
fsm = evstr(f+'(x+dx-dy)');
|
||||
fms = evstr(f+'(x-dx+dy)');
|
||||
fmm = evstr(f+'(x-dx-dy)');
|
||||
end
|
||||
D1 = (fss-fsm-fms+fmm)/(4*hx*hy);
|
||||
mask = ones(nf,1);
|
||||
j = 2;
|
||||
while (1)
|
||||
D2 = zeros(nf,j);
|
||||
hx = hx/con;
|
||||
hy = hy/con;
|
||||
dx = zeros(nx,1);
|
||||
dy = zeros(nx,1);
|
||||
dx(xi) = hx;
|
||||
dy(xj) = hy;
|
||||
if nargin > 3
|
||||
fss = evstr(f+'(x+dx+dy,varargin)');
|
||||
fsm = evstr(f+'(x+dx-dy,varargin)');
|
||||
fms = evstr(f+'(x-dx+dy,varargin)');
|
||||
fmm = evstr(f+'(x-dx-dy,varargin)');
|
||||
else
|
||||
fss = evstr(f+'(x+dx+dy)');
|
||||
fsm = evstr(f+'(x+dx-dy)');
|
||||
fms = evstr(f+'(x-dx+dy)');
|
||||
fmm = evstr(f+'(x-dx-dy)');
|
||||
end
|
||||
for fi = 1:nf
|
||||
if mask(fi)
|
||||
err2 = big;
|
||||
D2(fi,1) = (fss(fi)-fsm(fi)-fms(fi)+fmm(fi))/(4*hx*hy);
|
||||
fac = con2;
|
||||
for k = 2:j
|
||||
D2(fi,k) = D2(fi,k-1)+(D2(fi,k-1)-D1(fi,k-1))/(fac-1);
|
||||
fac = con2*fac;
|
||||
errt = max(abs(D2(fi,k)-D2(fi,k-1)),abs(D2(fi,k)-D1(fi,k-1)));
|
||||
if errt <= err2
|
||||
err2 = errt;
|
||||
deriv = D2(fi,k);
|
||||
end
|
||||
end
|
||||
err(fi) = abs(D2(fi,j)-D1(fi,j-1));
|
||||
if err(fi) < toler | err(fi) > safe*err2
|
||||
H(fi,(xi-1)*nx+xj) = deriv;
|
||||
H(fi,(xj-1)*nx+xi) = deriv;
|
||||
mask(fi) = 0;
|
||||
end
|
||||
end
|
||||
end
|
||||
if (mask == 0) then
|
||||
break
|
||||
end
|
||||
j = j+1;
|
||||
if j == maxit
|
||||
error('DIFFEXT didn''t converge. Try to increase gstep_ (default 0.01)')
|
||||
end
|
||||
D1 = D2;
|
||||
end
|
||||
[m_err,i] = max(err);
|
||||
if m_err > toler
|
||||
pause
|
||||
dyn_disp(D2)
|
||||
dyn_disp(err)
|
||||
dyn_disp([x dx dy])
|
||||
error('HESSEXT obtains an accuracy > 1e-12. Try to increase gstep_ (default 0.01)')
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
// 10/12/2001 MJ modified initial h
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
Fast hpfiltering:
|
||||
|
||||
[c,t]=hpfast(y,lambda)
|
||||
|
||||
Inputs:
|
||||
|
||||
y: vector (nx1) containing the series you wish to filter
|
||||
lambda: smoothing parameter (default: quarterly 1600)
|
||||
|
||||
Outputs:
|
||||
|
||||
c: Cyclical component (deviations from trend)
|
||||
t: Trend component
|
||||
|
||||
|
|
@ -0,0 +1,110 @@
|
|||
function [d,t]=hpfast(y,s)
|
||||
d=[];t=[];
|
||||
[nargout,nargin] = argn(0)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// Fast hpfiltering:
|
||||
//
|
||||
// [c,t]=hpfast(y,lambda)
|
||||
//
|
||||
// Inputs:
|
||||
//
|
||||
// y: vector (nx1) containing the series you wish to filter
|
||||
// lambda: smoothing parameter (default: quarterly 1600)
|
||||
//
|
||||
// Outputs:
|
||||
//
|
||||
// c: Cyclical component (deviations from trend)
|
||||
// t: Trend component
|
||||
//
|
||||
if nargin<2 then
|
||||
s = 1600;
|
||||
end
|
||||
iopt = 0;
|
||||
y = y(:);
|
||||
|
||||
//! unknown arg type, using mtlb_length
|
||||
n = mtlb_length(y);
|
||||
d1 = zeros(n,1);
|
||||
t = zeros(n,1);
|
||||
|
||||
v = zeros(3*n,6);
|
||||
if iopt~=1 then
|
||||
ss = s;
|
||||
nn = n;
|
||||
v11 = 1;
|
||||
v22 = 1;
|
||||
v12 = 0;
|
||||
for i = 3:n
|
||||
x = v11;
|
||||
z = v12;
|
||||
v11 = 1/s+4*(x-z)+v22;
|
||||
v12 = 2*x-z;
|
||||
v22 = x;
|
||||
det1 = v11*v22-v12*v12;
|
||||
v(i,1) = v22/det1;
|
||||
v(i,3) = v11/det1;
|
||||
v(i,2) = -v12/det1;
|
||||
x = v11+1;
|
||||
z = v11;
|
||||
v11 = v11-v11*v11/x;
|
||||
v22 = v22-v12*v12/x;
|
||||
v12 = v12-z*v12/x;
|
||||
end
|
||||
end
|
||||
//
|
||||
// Forward pass
|
||||
//
|
||||
m1 = y(2);
|
||||
m2 = y(1);
|
||||
for i = 3:n
|
||||
x = m1;
|
||||
m1 = 2*m1-m2;
|
||||
m2 = x;
|
||||
t(i-1,1) = v(i,1)*m1+v(i,2)*m2;
|
||||
d1(i-1,1) = v(i,2)*m1+v(i,3)*m2;
|
||||
det1 = v(i,1)*v(i,3)-v(i,2)*v(i,2);
|
||||
v11 = v(i,3)/det1;
|
||||
v12 = -v(i,2)/det1;
|
||||
z = (y(i)-m1)/(v11+1);
|
||||
m1 = m1+v11*z;
|
||||
m2 = m2+v12*z;
|
||||
end
|
||||
t(n,1) = m1;
|
||||
t(n-1,1) = m2;
|
||||
//
|
||||
// Backward pass
|
||||
//
|
||||
m1 = y(n-1);
|
||||
m2 = y(n);
|
||||
for i = n-2:-1:1
|
||||
i1 = i+1;
|
||||
ib = n-i+1;
|
||||
x = m1;
|
||||
m1 = 2*m1-m2;
|
||||
m2 = x;
|
||||
//
|
||||
// Combine info for y(.lt.i) with info for y(.ge.i)
|
||||
//
|
||||
if i>2 then
|
||||
e1 = v(ib,3)*m2+v(ib,2)*m1+t(i,1);
|
||||
e2 = v(ib,2)*m2+v(ib,1)*m1+d1(i,1);
|
||||
b11 = v(ib,3)+v(i1,1);
|
||||
b12 = v(ib,2)+v(i1,2);
|
||||
b22 = v(ib,1)+v(i1,3);
|
||||
det1 = b11*b22-b12*b12;
|
||||
t(i,1) = (-b12*e1+b11*e2)/det1;
|
||||
end
|
||||
//
|
||||
// End of combining
|
||||
//
|
||||
det1 = v(ib,1)*v(ib,3)-v(ib,2)*v(ib,2);
|
||||
v11 = v(ib,3)/det1;
|
||||
v12 = -v(ib,2)/det1;
|
||||
z = (y(i)-m1)/(v11+1);
|
||||
m1 = m1+v11*z;
|
||||
m2 = m2+v12*z;
|
||||
end
|
||||
t(1,1) = m1;
|
||||
t(2,1) = m2;
|
||||
d = y-t;
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
function [a]=indnv(x,y)
|
||||
a=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
%v=size(x)
|
||||
a = zeros(%v(1),%v(2));
|
||||
|
||||
for i = 1:size(x,1)
|
||||
j = find(x(i)==y)';
|
||||
if j==[] then
|
||||
a(i) = %nan;
|
||||
else
|
||||
a(i) = j(1);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
function [y_]=irf(dr, varexo, ssize, long_, drop_, replic, iorder)
|
||||
|
||||
if length(varexo) then
|
||||
shock_var = grep_exact(lgx_,varexo);
|
||||
else
|
||||
shock_var = 1;
|
||||
end
|
||||
|
||||
y_ = 0;
|
||||
for j = 1: replic
|
||||
rand('seed',j,'n');
|
||||
ex1_ = rand(long_+drop_+xkmin_,exo_nbr,'n')*chol(Sigma_e_);
|
||||
ex2_ = ex1_;
|
||||
ex2_(drop_+1,shock_var) = ex2_(drop_+1,shock_var)+ssize;
|
||||
y1_ = simult_(dr('ys'),dr,ex1_,iorder,long_+drop_);
|
||||
y2_ = simult_(dr('ys'),dr,ex2_,iorder,long_+drop_);
|
||||
y_ = y_+(y2_(:,ykmin_+drop_:$)-y1_(:,ykmin_+drop_:$));
|
||||
end
|
||||
y_ = y_/replic;
|
||||
dsample(1,long_-drop_);
|
||||
|
||||
dyn2vec();
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
function []=jacob(func,z)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
global('d1_','jacobia_','gstep_');
|
||||
|
||||
z = matrix(z,max(size(z)),1);
|
||||
d1_ = evstr(func+'(z)');
|
||||
nz = size(z,1);
|
||||
jacobia_ = zeros(size(d1_,1),nz);
|
||||
dh = max(abs(z),gstep_*ones(nz,1))*(%eps^(1/3));
|
||||
xdh = z;
|
||||
for j = 1:nz
|
||||
xdh(j) = xdh(j)+dh(j);
|
||||
h = xdh(j)-z(j);
|
||||
jacobia_(:,j) = (evstr(func+'(xdh)')-d1_)/h;
|
||||
xdh(j) = z(j);
|
||||
end
|
||||
|
||||
// 8/27/2000 modified with function name MJ
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
computes second order partial derivatives
|
||||
uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884
|
||||
|
|
@ -0,0 +1,74 @@
|
|||
function [jacob2_]=jacob2(func,x,varargin)
|
||||
jacob2_=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// computes second order partial derivatives
|
||||
// uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884
|
||||
|
||||
n = size(x,1);
|
||||
nargin = length(varargin);
|
||||
|
||||
//h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
|
||||
h1 = max(abs(x),sqrt(gstep_)*ones(n,1))*(%eps^(1/6));
|
||||
h_1 = h1;
|
||||
xh1 = x+h1;
|
||||
h1 = xh1-x;
|
||||
xh1 = x-h_1;
|
||||
h_1 = x-xh1;
|
||||
xh1 = x;
|
||||
|
||||
if nargin > 2 then
|
||||
f0 = evstr(func+'(x,varargin)');
|
||||
else
|
||||
f0 = evstr(func+'(x)');
|
||||
end
|
||||
f1 = zeros(size(f0,1),n);
|
||||
f_1 = f1;
|
||||
|
||||
for i = 1:n
|
||||
xh1(i) = x(i)+h1(i)
|
||||
if nargin > 2
|
||||
f1(:,i) = evstr(func+'(xh1,varargin)');
|
||||
else
|
||||
f1(:,i) = evstr(func+'(xh1)');
|
||||
end
|
||||
xh1(i) = x(i)-h_1(i);
|
||||
if nargin > 2
|
||||
f_1(:,i) = evstr(func+'(xh1,varargin)');
|
||||
else
|
||||
f_1(:,i) = evstr(func+'(xh1)');
|
||||
end
|
||||
xh1(i) = x(i);
|
||||
i = i+1;
|
||||
end
|
||||
xh_1 = xh1;
|
||||
|
||||
jacob2_ = sparse([],[],[size(f0,1) n*n]);
|
||||
for i = 1:n
|
||||
if i>1 then
|
||||
k = i:n:n*(i-1);
|
||||
jacob2_(:,(i-1)*n+1:(i-1)*n+i-1) = jacob2_(:,k);
|
||||
end
|
||||
jacob2_(:,(i-1)*n+i) = sparse((f1(:,i)+f_1(:,i)-2*f0) ./ (h1(i)*h_1(i)));
|
||||
temp = f1+f_1-f0*ones(1,n);
|
||||
for j = i+1:n
|
||||
xh1(i) = x(i)+h1(i);
|
||||
xh1(j) = x(j)+h_1(j);
|
||||
xh_1(i) = x(i)-h1(i);
|
||||
xh_1(j) = x(j)-h_1(j);
|
||||
if nargin > 2 then
|
||||
jacob2_(:,(i-1)*n+j) = sparse(-(-evstr(func+'(xh1,varargin)')-evstr(func+'(xh_1,varargin)')+temp(:,i)+temp(:,j)) ./ (2*h1(i)*h_1(j)));
|
||||
else
|
||||
jacob2_(:,(i-1)*n+j) = sparse(-(-evstr(func+'(xh1)')-evstr(func+'(xh_1)')+temp(:,i)+temp(:,j)) ./ (2*h1(i)*h_1(j)));
|
||||
end
|
||||
xh1(i) = x(i);
|
||||
xh1(j) = x(j);
|
||||
xh_1(i) = x(i);
|
||||
xh_1(j) = x(j);
|
||||
j = j+1;
|
||||
end
|
||||
i = i+1;
|
||||
end
|
||||
|
||||
// 10/21/02 MJ used the 7 points formula
|
||||
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
symmetric formula to compute the Jacobian
|
||||
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
function [jacobia_]=jacob_a(func,z,varargin)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// symmetric formula to compute the Jacobian
|
||||
//
|
||||
|
||||
global d1_
|
||||
nargin = length(varargin);
|
||||
|
||||
if nargin > 2 then
|
||||
d1_ = evstr(func+'(z,varargin)');
|
||||
else
|
||||
d1_ = evstr(func+'(z)');
|
||||
end
|
||||
nz = size(z,1);
|
||||
jacobia_ = zeros(size(d1_,1),nz);
|
||||
dh = max(abs(z),gstep_*ones(nz,1))*(%eps^(1/3));
|
||||
xdh1 = z;
|
||||
xdh2 = z;
|
||||
for j = 1:nz
|
||||
xdh1(j) = z(j)-dh(j);
|
||||
xdh2(j) = z(j)+dh(j);
|
||||
h = xdh2(j)-xdh1(j);
|
||||
if nargin > 2 then
|
||||
jacobia_(:,j) = (evstr(func+'(xdh2,varargin)')-evstr(func+'(xdh1,varargin)'))/h;
|
||||
else
|
||||
jacobia_(:,j) = (evstr(func+'(xdh2)')-evstr(func+'(xdh1)'))/h;
|
||||
end
|
||||
xdh1(j) = z(j);
|
||||
xdh2(j) = z(j);
|
||||
end
|
||||
|
||||
// 10/21/02 MJ creation
|
||||
Binary file not shown.
|
|
@ -0,0 +1,3 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
LINEAR : LINEAR (x,'filename')
|
||||
|
|
@ -0,0 +1,104 @@
|
|||
function []=linear(x)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// LINEAR : LINEAR (x,'filename')
|
||||
|
||||
|
||||
x = matrix(x,max(size(x)),1);
|
||||
xx = x;
|
||||
nn = size(iy_,2);
|
||||
nv = max(iy_(size(iy_,1),:));
|
||||
indic = zeros(0,0);
|
||||
|
||||
if size(x,1)==nv then
|
||||
r = matrix(find(matrix(iy_'>0,(ykmin_+ykmax_+1)*nn,1)),1,-1);
|
||||
n = abs(lgy_);
|
||||
n = ones(ykmin_+ykmax_+1,1).*.n;
|
||||
n = ascii(n(r,:));
|
||||
m = ((-ykmin_:ykmax_)')*ones(1,nn);
|
||||
m = matrix(m',(ykmin_+ykmax_+1)*nn,1);
|
||||
m = m(r);
|
||||
elseif size(x,1)==nn then
|
||||
l = iy_>0;
|
||||
l = bool2s(l) .* (ones(size(iy_,1),1).*.(1:nn));
|
||||
|
||||
//!! Unknown function nonzeros ,the original calling sequence is used
|
||||
i = nonzeros(matrix(iy_,size(iy_,1)*nn,1));
|
||||
|
||||
//!! Unknown function nonzeros ,the original calling sequence is used
|
||||
j = nonzeros(matrix(l,size(iy_,1)*nn,1));
|
||||
s = ones(1,max(iy_(size(iy_,1),:)));
|
||||
indic = sparse([i(:),j(:)],s,[nv,nn]);
|
||||
x = indic*x;
|
||||
n = lgy_;
|
||||
m = zeros(nn,1);
|
||||
else
|
||||
error('Wrong number of arguments in LINEAR.');
|
||||
end
|
||||
|
||||
jacob(x,'ff_');
|
||||
|
||||
if ~(indic==[]) then
|
||||
jacobia_ = jacobia_*indic;
|
||||
clear('indic');
|
||||
end
|
||||
|
||||
mtlb_fprintf(1,'Periods : ');
|
||||
mtlb_fprintf(1,'%4g \n',iter_);
|
||||
mtlb_fprintf(1,'Endogenous variables : ');
|
||||
mtlb_fprintf(1,'\n');
|
||||
mtlb_fprintf(1,lgy_);
|
||||
mtlb_fprintf(1,'\n');
|
||||
mtlb_fprintf(1,'Exogenous variables : ');
|
||||
mtlb_fprintf(1,'\n');
|
||||
mtlb_fprintf(1,lgx_);
|
||||
mtlb_fprintf(1,'\n');
|
||||
mtlb_fprintf(1,'Linearization around :');
|
||||
mtlb_fprintf(1,'\n');
|
||||
|
||||
for i = 1:size(n,1)
|
||||
mtlb_fprintf(1,n(i,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(i));
|
||||
mtlb_fprintf(1,' = %15.6f \n',xx(i));
|
||||
end
|
||||
|
||||
mtlb_fprintf(1,'\n');
|
||||
|
||||
for i = 1:size(jacobia_,1)
|
||||
for j = 1:size(jacobia_,2)
|
||||
if jacobia_(i,j)~=0 then
|
||||
if jacobia_(i,j)==1 then
|
||||
if j==1 then
|
||||
mtlb_fprintf(1,n(j,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(j));
|
||||
else
|
||||
mtlb_fprintf(1,' + ');
|
||||
mtlb_fprintf(1,n(j,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(j));
|
||||
end
|
||||
elseif jacobia_(i,j)==(-1) then
|
||||
mtlb_fprintf(1,' - ');
|
||||
mtlb_fprintf(1,n(j,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(j));
|
||||
elseif jacobia_(i,j)>0 then
|
||||
if j>1 then
|
||||
mtlb_fprintf(1,' + ');
|
||||
end
|
||||
mtlb_fprintf(1,'%15.6g',jacobia_(i,j));
|
||||
mtlb_fprintf(1,'*');
|
||||
mtlb_fprintf(1,n(j,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(j));
|
||||
else
|
||||
mtlb_fprintf(1,'%15.6g',jacobia_(i,j));
|
||||
mtlb_fprintf(1,'*');
|
||||
mtlb_fprintf(1,n(j,:));
|
||||
mtlb_fprintf(1,'(%1g)',m(j));
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
mtlb_fprintf(1,'\n');
|
||||
end
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
function [x,f,fvec,%check]=lnsrch(xold,fold,g,p,stpmax,func,varargin)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
alf = .0001;
|
||||
tolx = .000000000037;
|
||||
alam = 1;
|
||||
|
||||
nn = size(xold,1);
|
||||
|
||||
summ = sqrt(sum(p .* p));
|
||||
if abs(summ) > abs(stpmax) then
|
||||
p = p .* stpmax/summ;
|
||||
end
|
||||
|
||||
slope = g'*p;
|
||||
|
||||
test = max(abs(p)' ./ max([abs(xold)';ones(1,nn)],'r'));
|
||||
alamin = tolx/test;
|
||||
|
||||
if alamin>.1 then
|
||||
alamin = .1;
|
||||
end
|
||||
|
||||
while 1 then
|
||||
if alam<alamin then
|
||||
%check = 1;
|
||||
break
|
||||
;
|
||||
|
||||
end
|
||||
|
||||
x = xold+alam*p;
|
||||
|
||||
if argn(2) > 6
|
||||
fvec = evstr(func+'(x,varargin)');
|
||||
else
|
||||
fvec = evstr(func+'(x)');
|
||||
end
|
||||
|
||||
f = real(.5*fvec'*fvec);
|
||||
if or(isnan(fvec)) then
|
||||
alam = alam/2;
|
||||
alam2 = alam;
|
||||
f2 = f;
|
||||
fold2 = fold;
|
||||
else
|
||||
if abs(f) <= abs(fold+alf*alam*slope) then
|
||||
%check = 0;
|
||||
break;
|
||||
else
|
||||
if alam==1 then
|
||||
tmplam = -slope/(2*(f-fold-slope));
|
||||
else
|
||||
rhs1 = f-fold-alam*slope;
|
||||
rhs2 = f2-fold2-alam2*slope;
|
||||
a = (rhs1/(alam^2)-rhs2/(alam2^2))/(alam-alam2);
|
||||
b = (-alam2*rhs1/(alam^2)+alam*rhs2/(alam2^2))/(alam-alam2);
|
||||
if a==0 then
|
||||
tmplam = -slope/(2*b);
|
||||
else
|
||||
disc = b^2-3*a*slope;
|
||||
|
||||
if abs(disc)<0 then
|
||||
error('Roundoff problem in nlsearch');
|
||||
else
|
||||
tmplam = (-b+sqrt(disc))/(3*a);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if abs(tmplam) >.5*abs(alam) then
|
||||
tmplam = .5*alam;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
alam2 = alam;
|
||||
f2 = f;
|
||||
fold2 = fold;
|
||||
alam = max([abs(tmplam);.1*abs(alam)]);
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
// 01/14/01 MJ lnsearch is now a separate function
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
MCOMPARE : MCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
|
||||
This optional command plots the relative differences between
|
||||
two different simulations for a list of variables. One plot
|
||||
is drawn for each variable. The trajectories must have been
|
||||
previously saved by the instruction DYNASAVE. The simulation
|
||||
in file1 serves as the base simulation and the ploted quantity
|
||||
is equal to the difference between the two simulation reported
|
||||
to the first one. If, for a given variable, zero is one of the
|
||||
value of the base simulation, the absolute difference is ploted
|
||||
instead of the relative one.
|
||||
|
|
@ -0,0 +1,58 @@
|
|||
function []=mcompare(s1,s2)
|
||||
// MCOMPARE : MCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// This optional command plots the relative differences between
|
||||
// two different simulations for a list of variables. One plot
|
||||
// is drawn for each variable. The trajectories must have been
|
||||
// previously saved by the instruction DYNASAVE. The simulation
|
||||
// in file1 serves as the base simulation and the ploted quantity
|
||||
// is equal to the difference between the two simulation reported
|
||||
// to the first one. If, for a given variable, zero is one of the
|
||||
// value of the base simulation, the absolute difference is ploted
|
||||
// instead of the relative one.
|
||||
|
||||
global('dsmpl_','iter_');
|
||||
global('nvx','nvy','x','y','lag1');
|
||||
|
||||
ftest(s1,s2);
|
||||
|
||||
ix = (1-lag1(1):size(x,2)-lag1(1))';
|
||||
i = (lag1(1):size(ix,1)-lag1(2)+1)';
|
||||
|
||||
if size(dsmpl_,1)==1 then
|
||||
error('DSAMPLE not specified.');
|
||||
end
|
||||
|
||||
if dsmpl_(3)>0 then
|
||||
if dsmpl_(3)==2 then
|
||||
if dsmpl_(1)<0|dsmpl_(2)>(size(x,2)-lag1(2)) then
|
||||
error('Wrong sample.');
|
||||
end
|
||||
i = (dsmpl_(1)+lag1(1):dsmpl_(2)+lag1(1))';
|
||||
elseif dsmpl_(3)==1 then
|
||||
if dsmpl_(1)>(size(x,2)-lag1(2)) then
|
||||
error('Wrong sample.');
|
||||
end
|
||||
i = (lag1(1):dsmpl_(1)+lag1(1))';
|
||||
end
|
||||
end
|
||||
|
||||
for k = 1:size(x,1)
|
||||
xset('window',max(winsid()+1))
|
||||
max(winsid());
|
||||
x1 = x(k,i);
|
||||
y1 = y(k,i);
|
||||
|
||||
//! unknown arg type, using mtlb_length
|
||||
if nnz(x1)<mtlb_length(x1) then
|
||||
mtlb_plot(ix(i),y1-x1);
|
||||
else
|
||||
mtlb_plot(ix(i),(y1-x1) ./ x1);
|
||||
end
|
||||
xtitle(' ','Periods',' ');
|
||||
xtitle('Variable '+s2(k),' ',' ');
|
||||
end
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1,139 @@
|
|||
dynare
|
||||
disp_dr
|
||||
rplot
|
||||
simul
|
||||
check
|
||||
dyn_disp
|
||||
simul
|
||||
dyn_disp
|
||||
disp_th_moments
|
||||
stoch_simul
|
||||
irf
|
||||
disp_th_moments
|
||||
stoch_simul
|
||||
irf
|
||||
steady
|
||||
disp_moments
|
||||
disp_dr
|
||||
d_corr
|
||||
check
|
||||
dyn_mprintf
|
||||
simult
|
||||
ffill
|
||||
simult_
|
||||
bksup
|
||||
simk
|
||||
bksup1
|
||||
sim1
|
||||
bksupk
|
||||
shift
|
||||
brm
|
||||
selif
|
||||
bseastr
|
||||
sci_th_autocovariances
|
||||
uplib
|
||||
sci_shift
|
||||
linear
|
||||
sci_jacob_a
|
||||
dcompare
|
||||
sci_jacob2
|
||||
diffext
|
||||
sci_disp_th_moments
|
||||
rplot
|
||||
translatepaths
|
||||
th_autocovariances
|
||||
rfrot
|
||||
fff
|
||||
resol
|
||||
str2mat
|
||||
resid
|
||||
dr1
|
||||
reshapel
|
||||
dsample
|
||||
nonzeros
|
||||
dtranslatepaths
|
||||
mcompare
|
||||
dyn2vec
|
||||
lnsrch
|
||||
steady_
|
||||
fff_
|
||||
%s_2_s
|
||||
jacob_a
|
||||
ftest
|
||||
jacob2
|
||||
dynasave
|
||||
jacob
|
||||
solve
|
||||
dynatype
|
||||
equiv
|
||||
indnv
|
||||
fx_
|
||||
hpfast
|
||||
ff1_
|
||||
hessext
|
||||
ff2_
|
||||
grep_exact
|
||||
ff2a_
|
||||
gcompare
|
||||
ff_
|
||||
steady
|
||||
disp_moments
|
||||
d_corr
|
||||
dyn_mprintf
|
||||
simult
|
||||
fff_
|
||||
ffill
|
||||
fff
|
||||
ftest
|
||||
ff_
|
||||
fx_
|
||||
ff2a_
|
||||
gcompare
|
||||
ff2_
|
||||
grep_exact
|
||||
ff1_
|
||||
hessext
|
||||
equiv
|
||||
indnv
|
||||
dynatype
|
||||
solve
|
||||
dynasave
|
||||
jacob
|
||||
jacob2
|
||||
jacob_a
|
||||
linear
|
||||
steady_
|
||||
dyn2vec
|
||||
lnsrch
|
||||
dtranslatepaths
|
||||
mcompare
|
||||
dsample
|
||||
nonzeros
|
||||
dr1
|
||||
reshapel
|
||||
resid
|
||||
str2mat
|
||||
resol
|
||||
rfrot
|
||||
th_autocovariances
|
||||
simult_
|
||||
diffext
|
||||
sci_disp_th_moments
|
||||
dcompare
|
||||
sci_jacob2
|
||||
sci_jacob_a
|
||||
translatepaths
|
||||
sci_shift
|
||||
uplib
|
||||
bseastr
|
||||
sci_th_autocovariances
|
||||
brm
|
||||
selif
|
||||
bksupk
|
||||
shift
|
||||
bksup1
|
||||
sim1
|
||||
bksup
|
||||
simk
|
||||
%s_2_s
|
||||
hpfast
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
function y=nonzeros(x);
|
||||
x=x(:);
|
||||
k = find(x);
|
||||
y=x(k);
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
NOTES ABOUT DEVELOPMENT OF SCILAB VERSION OF DYNARE v. 4
|
||||
|
||||
|
||||
- we develop under scilab-4.0-b4-20060515
|
||||
- the first objective is to translate the estimation module
|
||||
- the two test cases are tests/fs2000/fs2000.mod (with numerical computation of the steadystate) and tests/fs2000/fs2000a.mod (with a *_steadystate.m model)
|
||||
- the dependencies to and from each Matlab routines (generated by m2html) is available in
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
RESHAPEL Change size.
|
||||
RESHAPEL(X,M,N) returns the M-by-N matrix whose elements
|
||||
are taken linewise from X. An error results if X does
|
||||
not have M*N elements.
|
||||
|
|
@ -0,0 +1,29 @@
|
|||
function [y]=reshapel(x,m,n)
|
||||
y=[];
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
//RESHAPEL Change size.
|
||||
// RESHAPEL(X,M,N) returns the M-by-N matrix whose elements
|
||||
// are taken linewise from X. An error results if X does
|
||||
// not have M*N elements.
|
||||
|
||||
[mm,nn] = size(x);
|
||||
|
||||
if mm*nn~=m*n then
|
||||
error('Matrix must have M*N elements.');
|
||||
end
|
||||
|
||||
%v = x';
|
||||
[%i,%i] = find(%v);%i=%i(:);%i=%i(:);
|
||||
%v = %v(:)
|
||||
if %i<>[] then s = %v(%i+size(%v,1)*(%i-1)) ;else s = [],end
|
||||
if size(i,2)~=1 then
|
||||
i = i';
|
||||
end
|
||||
k = (j-1)*nn+i;
|
||||
j = k-1-fix(k-1./n).*n+1;
|
||||
i = (k-j)/n+1;
|
||||
y = full(sparse([i(:),j(:)],s,[m,n]));
|
||||
|
||||
return
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,32 @@
|
|||
function []=resid()
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
global('iter_','valf_','ex_','y_','it_','exe_','ys_','iy_','ykmin_','ykmax_','endval_','z');
|
||||
|
||||
n = size(iy_,2);
|
||||
// if ~ valf_ | size(y_,2) ~= iter_+ykmin_+ykmax_
|
||||
if ~valf_ then
|
||||
if (size(ys_,1)==1)&(ys_==0) then
|
||||
ys_ = zeros(size(ys_,1),1);
|
||||
end
|
||||
y_ = ys_*ones(1,iter_+ykmin_+ykmax_);
|
||||
if endval_==1 then
|
||||
|
||||
//! mtlb(ys0_) can be replaced by ys0_() or ys0_ whether ys0_ is an m-file or not
|
||||
y_(:,1:ykmin_) = mtlb(ys0_)*ones(1,ykmin_);
|
||||
end
|
||||
end
|
||||
|
||||
i = iy_';
|
||||
iyr0 = find(i(:))';
|
||||
|
||||
y = y_(:);
|
||||
z = zeros(n,iter_);
|
||||
for it_ = ykmin_+1:iter_+ykmin_
|
||||
z(:,it_-ykmin_) = ff_(y(iyr0));
|
||||
iyr0 = iyr0+n;
|
||||
end
|
||||
|
||||
dyn_disp([(1:iter_)',z']);
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
function [dr_]=resol(ys,algo,%linear,iorder)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
|
||||
xlen = xkmax_+xkmin_+1;
|
||||
klen = ykmin_+ykmax_+1;
|
||||
iyv = iy_';
|
||||
iyv = iyv(:);
|
||||
iyr0 = find(iyv)';
|
||||
it_ = ykmin_+1;
|
||||
|
||||
|
||||
if exo_nbr==0 then
|
||||
exe_ = [];
|
||||
end
|
||||
|
||||
if ~(iy_(ykmin_+1,:)>0) then
|
||||
error('RESOL: Error in model specification: some variables don""t appear as current');
|
||||
end
|
||||
|
||||
if ykmax_==0 then
|
||||
error('Backward or static model: no point in using RESOL');
|
||||
end
|
||||
|
||||
if xlen>1 then
|
||||
error('RESOL: stochastic exogenous variables must appear only at the'+' current period. Use additional endogenous variables');
|
||||
end
|
||||
|
||||
dr_ = struct();
|
||||
// check if ys is steady state
|
||||
fff_name = fname_+'_fff';
|
||||
if max(abs(evstr(fff_name+'(ys)'))) > dynatol_ then
|
||||
[ys_temp,%check]=solve(fff_name,ys);
|
||||
dr_('ys') = ys_temp;
|
||||
if %check then
|
||||
error('RESOL: convergence problem in SOLVE');
|
||||
end
|
||||
else
|
||||
dr_('ys') = ys
|
||||
end
|
||||
|
||||
dr_ = dr1(iorder,dr_);
|
||||
|
||||
// 08/24/2001 MJ uses Schmitt-Grohe and Uribe (2001) constant correction
|
||||
// in dr.ghs2
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
apply fast Givens rotation to a pair of rows
|
||||
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
function [a]=rfrot(a,alph,bet,%type)
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// apply fast Givens rotation to a pair of rows
|
||||
//
|
||||
if %type==1 then
|
||||
tau = bet*a(1,:)+a(2,:);
|
||||
a(2,:) = a(1,:)+alph*a(2,:);
|
||||
a(1,:) = tau;
|
||||
else
|
||||
tau = a(1,:)+bet*a(2,:);
|
||||
a(2,:) = alph*a(1,:)+a(2,:);
|
||||
a(1,:) = tau;
|
||||
end
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
Copyright (C) 2001 Michel Juillard
|
||||
|
||||
RPLOT : RPLOT ( ['var1'; 'var2'; ...] , rplottype )
|
||||
This optionnal command creates the plot of the variable
|
||||
trajectory. By default, the entire simulation period is
|
||||
ploted. The instruction DSAMPLE permits to reduce the
|
||||
number of periods in the plot.
|
||||
|
|
@ -0,0 +1,86 @@
|
|||
function []=rplot(s1,rplottype)
|
||||
[nargout,nargin] = argn(0)
|
||||
// RPLOT : RPLOT ( ['var1'; 'var2'; ...] , rplottype )
|
||||
// Copyright (C) 2001 Michel Juillard
|
||||
//
|
||||
// This optionnal command creates the plot of the variable
|
||||
// trajectory. By default, the entire simulation period is
|
||||
// ploted. The instruction DSAMPLE permits to reduce the
|
||||
// number of periods in the plot.
|
||||
|
||||
global dsmpl_ y_
|
||||
|
||||
if nargin==1 then
|
||||
rplottype = 0;
|
||||
end
|
||||
|
||||
ix = (1-ykmin_:size(y_,2)-ykmin_)';
|
||||
|
||||
m = size(s1,1);
|
||||
k=[];
|
||||
for j=1:m
|
||||
k1 = grep_exact(lgy_,s1(j));
|
||||
if k1==[] then
|
||||
error('One of the variable specified does not exist');
|
||||
end
|
||||
k = [k ; k1];
|
||||
end
|
||||
|
||||
|
||||
y = y_(k,:)';
|
||||
if dsmpl_==0 then
|
||||
i = (ykmin_:size(y_,2))';
|
||||
else
|
||||
i = (dsmpl_(1):dsmpl_(2))';
|
||||
end
|
||||
|
||||
t = 'Plot of ';
|
||||
if rplottype==0 then
|
||||
xset('window',max(winsid()+1))
|
||||
xset('use color',0)
|
||||
leg = "";
|
||||
for j=1:m
|
||||
t = t+s1(j)+' ';
|
||||
if j > 1
|
||||
leg = leg+'@'
|
||||
end
|
||||
leg = leg+s1(j);
|
||||
end
|
||||
plot2d(ix(i),y(i,:),leg=leg);
|
||||
xtitle(t,'Periods',' ');
|
||||
elseif rplottype==1 then
|
||||
for j = 1:size(y,1)
|
||||
xset('window',max(winsid()+1))
|
||||
max(winsid());
|
||||
plot2d(ix(i),y(j,i));
|
||||
xtitle('Plot of '+s1(:,j),' ',' ');
|
||||
xtitle(' ','Periods',' ');
|
||||
end
|
||||
elseif rplottype==2 then
|
||||
xset('window',max(winsid()+1))
|
||||
max(winsid());
|
||||
nl = max(1,fix(size(y,1)/4));
|
||||
nc = ceil(size(y,1)/nl);
|
||||
for j = 1:size(y,1)
|
||||
%v2$2 = int((j-1)/nc);%v2$1 = j-1-nc*%v2$2
|
||||
xsetech([%v2$1/nc,%v2$2/nl,1/nc,1/nl]);
|
||||
mtlb_plot(ix(i),y(j,i));
|
||||
mtlb_hold('on');
|
||||
mtlb_plot(ix(i),ys_(j)*ones(1,size(i,1)),'w:');
|
||||
xtitle(' ','Periods',' ');
|
||||
xtitle(' ',' ',s1(:,j));
|
||||
xtitle('Plot of '+s1(:,j),' ',' ');
|
||||
end
|
||||
end
|
||||
|
||||
// 02/28/01 MJ replaced bseastr by MATLAB's strmatch
|
||||
// 06/19/01 MJ added 'exact' to strmatch calls
|
||||
// 09/25/01 MJ translated to Scilab
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
function [stk,txt,top]=sci_disp_th_moments()
|
||||
RHS=[]
|
||||
for k=1:rhs
|
||||
RHS=[stk(top)(1),RHS]
|
||||
top=top-1
|
||||
end
|
||||
top=top+1
|
||||
stk=list('disp_th_moments'+rhsargs(RHS),'0','?','?','?')
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
function [stk,txt,top]=sci_jacob2()
|
||||
RHS=[]
|
||||
for k=1:rhs
|
||||
RHS=[stk(top)(1),RHS]
|
||||
top=top-1
|
||||
end
|
||||
top=top+1
|
||||
stk=list('jacob2'+rhsargs(RHS),'0','?','?','?')
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue