v4: cleanup unused function
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@780 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
michel
parent
757a5f4af2
commit
ac5a09a233
|
|
@ -1,186 +0,0 @@
|
|||
function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,gend,data)
|
||||
% stephane.adjemian@cepremap.cnrs.fr [09-07-2004]
|
||||
%
|
||||
% Adapted from mj_optmumlik.m
|
||||
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_ xparam1_test
|
||||
|
||||
fval = [];
|
||||
ys = [];
|
||||
trend_coeff = [];
|
||||
xparam1_test = xparam1;
|
||||
cost_flag = 1;
|
||||
nobs = size(options_.varobs,1);
|
||||
%------------------------------------------------------------------------------
|
||||
% 1. Get the structural parameters & define penalties
|
||||
%------------------------------------------------------------------------------
|
||||
if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
|
||||
k = find(xparam1 < bayestopt_.lb);
|
||||
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return;
|
||||
end
|
||||
if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
|
||||
k = find(xparam1 > bayestopt_.ub);
|
||||
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return;
|
||||
end
|
||||
Q = M_.Sigma_e;
|
||||
for i=1:estim_params_.nvx
|
||||
k =estim_params_.var_exo(i,1);
|
||||
Q(k,k) = xparam1(i)*xparam1(i);
|
||||
end
|
||||
offset = estim_params_.nvx;
|
||||
if estim_params_.nvn
|
||||
H = zeros(nobs,nobs);
|
||||
for i=1:estim_params_.nvn
|
||||
k = estim_params_.var_endo(i,1);
|
||||
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
|
||||
end
|
||||
offset = offset+estim_params_.nvn;
|
||||
end
|
||||
if estim_params_.ncx
|
||||
for i=1:estim_params_.ncx
|
||||
k1 =estim_params_.corrx(i,1);
|
||||
k2 =estim_params_.corrx(i,2);
|
||||
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
|
||||
Q(k2,k1) = Q(k1,k2);
|
||||
end
|
||||
[CholQ,testQ] = chol(Q);
|
||||
if testQ %% The variance-covariance matrix of the structural innovations is not definite positive.
|
||||
%% We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||
a = diag(eig(Q));
|
||||
k = find(a < 0);
|
||||
if k > 0
|
||||
fval = bayestopt_.penalty+sum(-a(k));
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return
|
||||
end
|
||||
end
|
||||
offset = offset+estim_params_.ncx;
|
||||
end
|
||||
if estim_params_.ncn
|
||||
for i=1:estim_params_.ncn
|
||||
k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
|
||||
k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
|
||||
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
|
||||
H(k2,k1) = H(k1,k2);
|
||||
end
|
||||
[CholH,testH] = chol(H);
|
||||
if testH
|
||||
a = diag(eig(H));
|
||||
k = find(a < 0);
|
||||
if k > 0
|
||||
fval = bayestopt_.penalty+sum(-a(k));
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return
|
||||
end
|
||||
end
|
||||
offset = offset+estim_params_.ncn;
|
||||
end
|
||||
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
|
||||
% for i=1:estim_params_.np
|
||||
% M_.params(estim_params_.param_vals(i,1)) = xparam1(i+offset);
|
||||
%end
|
||||
M_.Sigma_e = Q;
|
||||
%------------------------------------------------------------------------------
|
||||
% 2. call model setup & reduction program
|
||||
%------------------------------------------------------------------------------
|
||||
[T,R,SteadyState,info] = dynare_resolve;
|
||||
if info(1) == 1 | info(1) == 2 | info(1) == 5
|
||||
fval = bayestopt_.penalty+1;
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return
|
||||
elseif info(1) == 3 | info(1) == 4 | info(1) == 20
|
||||
fval = bayestopt_.penalty+info(2)^2;
|
||||
llik=fval;
|
||||
cost_flag = 0;
|
||||
return
|
||||
end
|
||||
if options_.loglinear == 1
|
||||
constant = log(SteadyState(bayestopt_.mfys));
|
||||
else
|
||||
constant = SteadyState(bayestopt_.mfys);
|
||||
end
|
||||
if bayestopt_.with_trend == 1
|
||||
trend_coeff = zeros(nobs,1);
|
||||
for i=1:nobs
|
||||
trend_coeff(i) = evalin('base',bayestopt_.trend_coeff{i});
|
||||
end
|
||||
trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
|
||||
else
|
||||
trend = repmat(constant,1,gend);
|
||||
end
|
||||
start = options_.presample+1;
|
||||
np = size(T,1);
|
||||
mf = bayestopt_.mf;
|
||||
%------------------------------------------------------------------------------
|
||||
% 3. Initial condition of the Kalman filter
|
||||
%------------------------------------------------------------------------------
|
||||
if options_.lik_init == 1 % Kalman filter
|
||||
Pstar = lyapunov_symm(T,R*Q*transpose(R));
|
||||
Pinf = [];
|
||||
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
|
||||
Pstar = 10*eye(np);
|
||||
Pinf = [];
|
||||
elseif options_.lik_init == 3 % Diffuse Kalman filter
|
||||
Pstar = zeros(np,np);
|
||||
ivs = bayestopt_.i_T_var_stable;
|
||||
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R(ivs,:)*Q* ...
|
||||
transpose(R(ivs,:)));
|
||||
Pinf = bayestopt_.Pinf;
|
||||
% by M. Ratto
|
||||
RR=T(:,find(~ismember([1:np],ivs)));
|
||||
i=find(abs(RR)>1.e-10);
|
||||
R0=zeros(size(RR));
|
||||
R0(i)=sign(RR(i));
|
||||
Pinf=R0*R0';
|
||||
% by M. Ratto
|
||||
end
|
||||
%------------------------------------------------------------------------------
|
||||
% 4. Likelihood evaluation
|
||||
%------------------------------------------------------------------------------
|
||||
if estim_params_.nvn
|
||||
if options_.kalman_algo == 1
|
||||
[LIK, lik] = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start);
|
||||
if isinf(LIK) & ~estim_params_.ncn %% The univariate approach considered here doesn't
|
||||
%% apply when H has some off-diagonal elements.
|
||||
[LIK, lik] = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start);
|
||||
elseif isinf(LIK) & estim_params_.ncn
|
||||
[LIK, lik] = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start);
|
||||
end
|
||||
elseif options_.kalman_algo == 3
|
||||
if ~estim_params_.ncn %% The univariate approach considered here doesn't
|
||||
%% apply when H has some off-diagonal elements.
|
||||
[LIK, lik] = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start);
|
||||
else
|
||||
[LIK, lik] = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start);
|
||||
end
|
||||
end
|
||||
else
|
||||
if options_.kalman_algo == 1
|
||||
[LIK, lik] = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start);
|
||||
if isinf(LIK)
|
||||
[LIK, lik] = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start);
|
||||
end
|
||||
elseif options_.kalman_algo == 3
|
||||
[LIK, lik] = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start);
|
||||
end
|
||||
end
|
||||
if imag(LIK) ~= 0
|
||||
likelihood = bayestopt_.penalty;
|
||||
lik=ones(size(lik)).*bayestopt_.penalty;
|
||||
else
|
||||
likelihood = LIK;
|
||||
end
|
||||
% ------------------------------------------------------------------------------
|
||||
% Adds prior if necessary
|
||||
% ------------------------------------------------------------------------------
|
||||
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
|
||||
fval = (likelihood-lnprior);
|
||||
llik=[-lnprior; .5*lik(start:end)];
|
||||
|
|
@ -1,8 +1,8 @@
|
|||
function check_model()
|
||||
global iy_ ykmin_ ykmax_ xkmin_ xkmax_ iy_ exo_det_nbr
|
||||
global M_
|
||||
|
||||
xlen = xkmin_ + xkmax_ + 1;
|
||||
if ~ iy_(ykmin_+1,:) > 0
|
||||
xlen = M_.maximum_exo_lag+M_.maximum_exo_lead + 1;
|
||||
if ~ M_.lead_lag_incidence(M_.maximum_lag+1,:) > 0
|
||||
error ('RESOL: Error in model specification: some variables don"t appear as current') ;
|
||||
end
|
||||
|
||||
|
|
@ -11,7 +11,7 @@ if xlen > 1
|
|||
' current period. Use additional endogenous variables']) ;
|
||||
end
|
||||
|
||||
if (exo_det_nbr > 0) & (ykmin_ > 1 | ykmax_ > 1)
|
||||
if (M_.exo_det_nbr > 0) & (M_.maximum_lag > 1 | M_.maximum_lead > 1)
|
||||
error(['Exogenous deterministic variables are currently only allowed in' ...
|
||||
' models with leads and lags on only one period'])
|
||||
end
|
||||
|
|
|
|||
32
matlab/dr1.m
32
matlab/dr1.m
|
|
@ -45,7 +45,9 @@ z = z(iyr0) ;
|
|||
if options_.order == 1
|
||||
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,tempex);
|
||||
elseif options_.order == 2
|
||||
[junk,jacobia_,hessian] = feval([M_.fname '_dynamic'],z,tempex);
|
||||
[junk,jacobia_,hessian] = feval([M_.fname '_dynamic'],z,...
|
||||
[oo_.exo_simul ...
|
||||
oo_.exo_det_simul]);
|
||||
end
|
||||
|
||||
oo_.exo_simul = tempex ;
|
||||
|
|
@ -182,7 +184,7 @@ if any(isinf(a(:)))
|
|||
return
|
||||
end
|
||||
if M_.exo_nbr
|
||||
fu = aa(:,nz+1:end);
|
||||
fu = aa(:,nz+(1:M_.exo_nbr));
|
||||
end
|
||||
clear aa;
|
||||
|
||||
|
|
@ -317,15 +319,15 @@ if use_qzdiv
|
|||
end
|
||||
|
||||
%exogenous deterministic variables
|
||||
if M_.exo_det_nbr > 1
|
||||
if M_.exo_det_nbr > 0
|
||||
f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_lag+2:end,order_var))));
|
||||
f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_lag+1,order_var))));
|
||||
fudet = sparse(jacobia_(:,nz+M_.exo_nbr+1:end));
|
||||
M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nyf-nboth)]);
|
||||
M2 = M1*f1;
|
||||
dr.ghud = cell(M_.oo_.exo_simuldet_length,1);
|
||||
dr.ghud = cell(M_.exo_det_length,1);
|
||||
dr.ghud{1} = -M1*fudet;
|
||||
for i = 2:M_.oo_.exo_simuldet_length
|
||||
for i = 2:M_.exo_det_length
|
||||
dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:);
|
||||
end
|
||||
end
|
||||
|
|
@ -344,14 +346,14 @@ if M_.maximum_lag > 0
|
|||
end
|
||||
kk = kk';
|
||||
kk = find(kk(:));
|
||||
nk = size(kk,1)+M_.exo_nbr;
|
||||
nk = size(kk,1) + M_.exo_nbr + M_.exo_det_nbr;
|
||||
k1 = M_.lead_lag_incidence(:,order_var);
|
||||
k1 = k1';
|
||||
k1 = k1(:);
|
||||
k1 = k1(kk);
|
||||
k2 = find(k1);
|
||||
kk1(k1(k2)) = k2;
|
||||
kk1 = [kk1 length(k1)+1:length(k1)+M_.exo_nbr];
|
||||
kk1 = [kk1 length(k1)+1:length(k1)+M_.exo_nbr+M_.exo_det_nbr];
|
||||
kk = reshape([1:nk^2],nk,nk);
|
||||
kk1 = kk(kk1,kk1);
|
||||
%[junk,junk,hessian] = feval([M_.fname '_dynamic'],z, oo_.exo_steady_state);
|
||||
|
|
@ -384,8 +386,8 @@ for i=1:M_.maximum_lead
|
|||
kk = kk(:,k0);
|
||||
k0 = k1;
|
||||
end
|
||||
zx=[zx; zeros(M_.exo_nbr,np)];
|
||||
zu=[zu; eye(M_.exo_nbr)];
|
||||
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)];
|
||||
[n1,n2] = size(zx);
|
||||
if n1*n1*n2*n2 > 1e7
|
||||
rhs = zeros(M_.endo_nbr,n2*n2);
|
||||
|
|
@ -571,33 +573,33 @@ if M_.exo_det_nbr > 0
|
|||
hud = dr.ghud{1}(nstatic+1:nstatic+npred,:);
|
||||
zud=[zeros(np,M_.exo_det_nbr);dr.ghud{1};gx(:,1:npred)*hud;zeros(M_.exo_nbr,M_.exo_det_nbr);eye(M_.exo_det_nbr)];
|
||||
R1 = hessian*kron(zx,zud);
|
||||
dr.ghxud = cell(M_.oo_.exo_simuldet_length,1);
|
||||
dr.ghxud = cell(M_.exo_det_length,1);
|
||||
kf = [M_.endo_nbr-nyf+1:M_.endo_nbr];
|
||||
kp = nstatic+[1:npred];
|
||||
dr.ghxud{1} = -M1*(R1+f1*dr.ghxx(kf,:)*kron(dr.ghx(kp,:),dr.ghud{1}(kp,:)));
|
||||
Eud = eye(M_.exo_det_nbr);
|
||||
for i = 2:M_.oo_.exo_simuldet_length
|
||||
for i = 2:M_.exo_det_length
|
||||
hudi = dr.ghud{i}(kp,:);
|
||||
zudi=[zeros(np,M_.exo_det_nbr);dr.ghud{i};gx(:,1:npred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
|
||||
R2 = hessian*kron(zx,zudi);
|
||||
dr.ghxud{i} = -M2*(dr.ghxud{i-1}(kf,:)*kron(hx,Eud)+dr.ghxx(kf,:)*kron(dr.ghx(kp,:),dr.ghud{i}(kp,:)))-M1*R2;
|
||||
end
|
||||
R1 = hessian*kron(zu,zud);
|
||||
dr.ghudud = cell(M_.oo_.exo_simuldet_length,1);
|
||||
dr.ghudud = cell(M_.exo_det_length,1);
|
||||
kf = [M_.endo_nbr-nyf+1:M_.endo_nbr];
|
||||
|
||||
dr.ghuud{1} = -M1*(R1+f1*dr.ghxx(kf,:)*kron(dr.ghu(kp,:),dr.ghud{1}(kp,:)));
|
||||
Eud = eye(M_.exo_det_nbr);
|
||||
for i = 2:M_.oo_.exo_simuldet_length
|
||||
for i = 2:M_.exo_det_length
|
||||
hudi = dr.ghud{i}(kp,:);
|
||||
zudi=[zeros(np,M_.exo_det_nbr);dr.ghud{i};gx(:,1:npred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
|
||||
R2 = hessian*kron(zu,zudi);
|
||||
dr.ghuud{i} = -M2*dr.ghxud{i-1}(kf,:)*kron(hu,Eud)-M1*R2;
|
||||
end
|
||||
R1 = hessian*kron(zud,zud);
|
||||
dr.ghudud = cell(M_.oo_.exo_simuldet_length,M_.oo_.exo_simuldet_length);
|
||||
dr.ghudud = cell(M_.exo_det_length,M_.exo_det_length);
|
||||
dr.ghudud{1,1} = -M1*R1-M2*dr.ghxx(kf,:)*kron(hud,hud);
|
||||
for i = 2:M_.oo_.exo_simuldet_length
|
||||
for i = 2:M_.exo_det_length
|
||||
hudi = dr.ghud{i}(nstatic+1:nstatic+npred,:);
|
||||
zudi=[zeros(np,M_.exo_det_nbr);dr.ghud{i};gx(:,1:npred)*hudi+dr.ghud{i-1}(kf,:);zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
|
||||
R2 = hessian*kron(zudi,zudi);
|
||||
|
|
|
|||
Binary file not shown.
|
|
@ -1,511 +0,0 @@
|
|||
function hessian2(xparam1,gend,data)
|
||||
% source Schorfheide
|
||||
%/* filename: lbdhess.g
|
||||
%** description: The program computes the hessianfit at the posterior mode
|
||||
%** created: 05/05/00
|
||||
%/********************************************************
|
||||
%/* Compute hessianfit, element by element, fine tune with dxscale
|
||||
% Compute hessianfit for two step-seize (dx) and take average to prevent singularity
|
||||
%/********************************************************
|
||||
npara = length(xparam1);
|
||||
para = xparam1;
|
||||
ndx = 1;%6
|
||||
%dx = exp(-seqa(6,2,ndx));
|
||||
%dx = exp([-6:-2:-16]);
|
||||
dx = exp([-10]);
|
||||
hessianfit = zeros( npara, npara );
|
||||
gradx = zeros(ndx,1);
|
||||
grady = zeros(ndx,1);
|
||||
gradxy = zeros(ndx,1);
|
||||
hessdiag = zeros(ndx,1);
|
||||
dxscale = ones(npara,1);
|
||||
% dxscale(5,1)=10;
|
||||
% dxscale(13,1)=10;
|
||||
% dxscale(17,1)=10;
|
||||
% dxscale(8,1)=.10;
|
||||
% dxscale(11,1)=.10;
|
||||
% dxscale(31,1)=.10;
|
||||
|
||||
|
||||
%/* Compute Diagonal elements first
|
||||
%*/
|
||||
seli = 1;
|
||||
fx = mj_optmumlik(para,gend,data,1);
|
||||
%do until seli > npara;
|
||||
while seli <= npara;
|
||||
% locate 1,1;
|
||||
% "hessianfit Element (" seli seli ")";
|
||||
i=1;
|
||||
while i <= ndx;
|
||||
paradx = para;
|
||||
parady = para;
|
||||
paradx(seli) = paradx(seli) + dx(i)*dxscale(seli);
|
||||
parady(seli) = parady(seli) - dx(i)*dxscale(seli);
|
||||
paradxdy = paradx;
|
||||
paradxdy(seli) = paradxdy(seli) - dx(i)*dxscale(seli);
|
||||
% fx = optmumlik20(para,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdx = mj_optmumlik(paradx,gend,data,1);
|
||||
% fdx = optmumlik20(paradx,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdy = mj_optmumlik(parady,gend,data,1);
|
||||
% fdy = optmumlik20(parady,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
% fdxdy = optmumlik20(paradxdy,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdxdy = fx;
|
||||
gradx(i) = -( fx - fdx )/ (dx(i)*dxscale(seli));
|
||||
grady(i) = ( fx - fdy )/ (dx(i)*dxscale(seli));
|
||||
gradxy(i) = -(fx -fdxdy)/ sqrt( (dx(i)*dxscale(seli))^2 + (dx(i)*dxscale(seli))^2 );
|
||||
hessdiag(i) = -( 2*fx - fdx - fdy)/(dx(i)*dxscale(seli))^2;
|
||||
hessdiag(i) = -( fx - fdx - fdy + fdxdy )/(dx(i)*dx(i)*dxscale(seli)*dxscale(seli));
|
||||
i = i+1;
|
||||
end;
|
||||
% "Values";
|
||||
% -hessdiag';
|
||||
hessianfit(seli,seli) = -1*(hessdiag(1));
|
||||
% hessianfit(seli,seli) = -0.5*(hessdiag(3)+hessdiag(4));
|
||||
% locate 6,1;
|
||||
% "Value Used:" hessianfit[seli,seli];
|
||||
seli = seli+1
|
||||
end;
|
||||
|
||||
diag(hessianfit)
|
||||
|
||||
%/* Now compute off-diagonal elements
|
||||
%** Make sure that correlations are between -1 and 1
|
||||
%** errorij contains the index of elements that are invalid
|
||||
%*/
|
||||
errorij = [ 0 0 0];
|
||||
|
||||
seli = 1;
|
||||
for seli = 1:npara;
|
||||
% selj = seli+1;
|
||||
for selj =seli+1:npara;
|
||||
disp([seli selj]);
|
||||
% locate 1,1;
|
||||
% "hessianfit Element (" seli selj ")";
|
||||
i=1;
|
||||
while i <= ndx;
|
||||
paradx = para;
|
||||
parady = para;
|
||||
paradx(seli) = paradx(seli) + dx(i)*dxscale(seli);
|
||||
parady(selj) = parady(selj) - dx(i)*dxscale(selj);
|
||||
paradxdy = paradx;
|
||||
paradxdy(selj) = paradxdy(selj) - dx(i)*dxscale(selj);
|
||||
% fx = optmumlik20(para,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdx = mj_optmumlik(paradx,gend,data,1);
|
||||
% fdx = optmumlik20(paradx,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdy = mj_optmumlik(parady,gend,data,1);
|
||||
% fdy = optmumlik20(parady,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdy = mj_optmumlik(paradxdy,gend,data,1);
|
||||
% fdxdy = optmumlik20(paradxdy,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
gradx(i) = -( fx - fdx )/ (dx(i)*dxscale(seli));
|
||||
grady(i) = ( fx - fdy )/ (dx(i)*dxscale(selj));
|
||||
gradxy(i) = -(fx -fdxdy)/ sqrt( (dx(i)*dxscale(selj))^2 + (dx(i)*dxscale(seli))^2 );
|
||||
hessdiag(i) = -( 2*fx - fdx - fdy)/(dx(i)*dxscale(seli))^2;
|
||||
hessdiag(i) = -( fx - fdx - fdy + fdxdy )/(dx(i)*dx(i)*dxscale(seli)*dxscale(selj));
|
||||
i = i+1;
|
||||
end;
|
||||
% "Values";
|
||||
% -hessdiag';
|
||||
|
||||
% hessianfit(seli,selj) = -0.5*(hessdiag(3)+hessdiag(4));
|
||||
hessianfit(seli,selj) = -1*(hessdiag(1));
|
||||
|
||||
if ( hessianfit(seli,seli) == 0 ) | (hessianfit(selj,selj) == 0);
|
||||
corrij = 0;
|
||||
else;
|
||||
corrij = hessianfit(seli,selj)/sqrt(hessianfit(seli,seli)*hessianfit(selj,selj));
|
||||
end;
|
||||
|
||||
if (corrij < -1) | (corrij > 1);
|
||||
hessianfit(seli,selj)=0;
|
||||
errorij = [ errorij [seli selj corrij] ];
|
||||
end;
|
||||
hessianfit(selj,seli) = hessianfit(seli,selj);
|
||||
|
||||
% locate 6,1;
|
||||
% "Value Used: " hessianfit[seli,selj];
|
||||
% "Correlation:" corrij;
|
||||
% "Number of Errors:" rows(errorij)-1;
|
||||
% selj=selj+1;
|
||||
end;
|
||||
% seli = seli+1;
|
||||
end;
|
||||
|
||||
%cls;
|
||||
disp('Errors')
|
||||
disp(errorij);
|
||||
|
||||
|
||||
%/*******************************************************************************
|
||||
|
||||
bbbb=xparam1;
|
||||
% func =fval;
|
||||
% grad=grad;
|
||||
% retcode=exitflag
|
||||
opfhessfit = (-hessianfit);
|
||||
invhess=inv(opfhessfit);
|
||||
stdh=sqrt(diag(invhess));
|
||||
pr =length(xparam1);
|
||||
tstath=zeros(pr,1);
|
||||
i = 1; while i <= pr ; %do until i>pr;
|
||||
tstath(i)=abs(bbbb(i))/stdh(i);
|
||||
i=i+1; end ; %endo;
|
||||
%tstath
|
||||
% print "t-stats. from the Hessian";
|
||||
disp('print "t-stats. from the Hessian" ') ;
|
||||
disp([xparam1 stdh tstath]);
|
||||
|
||||
bbbb=xparam1;
|
||||
% func =fval;
|
||||
% grad=grad;
|
||||
% retcode=exitflag
|
||||
% opfhessfit = (-hessianfit);
|
||||
% invhess=inv(opfhessfit*.5+opfhess*.5);
|
||||
% stdh=sqrt(diag(invhess));
|
||||
% pr =length(xparam1);
|
||||
% tstath=zeros(pr,1);
|
||||
% i = 1; while i <= pr ; %do until i>pr;
|
||||
% tstath(i)=abs(bbbb(i))/stdh(i);
|
||||
% i=i+1; end ; %endo;
|
||||
%tstath
|
||||
% print "t-stats. from the Hessian";
|
||||
% disp('print "t-stats. from the Hessian" ') ;
|
||||
% disp([xparam1 stdh tstath]);
|
||||
|
||||
%opfhessfit = -hessianfit*.5+opfhess*.5;
|
||||
hessian=opfhessfit;
|
||||
|
||||
|
||||
|
||||
return;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
npara = length(xparam1);
|
||||
para = xparam1;
|
||||
ndx = 1;%6
|
||||
%dx = exp(-seqa(6,2,ndx));
|
||||
%dx = exp([-6:-2:-16]);
|
||||
dx = exp([-10]);
|
||||
hessianfit = zeros( npara, npara );
|
||||
gradx = zeros(ndx,1);
|
||||
grady = zeros(ndx,1);
|
||||
gradxy = zeros(ndx,1);
|
||||
hessdiag = zeros(ndx,1);
|
||||
dxscale = ones(npara,1);
|
||||
|
||||
|
||||
%/* Compute Diagonal elements first
|
||||
%*/
|
||||
seli = 1;
|
||||
%do until seli > npara;
|
||||
while seli <= npara;
|
||||
% locate 1,1;
|
||||
% "hessianfit Element (" seli seli ")";
|
||||
i=1;
|
||||
while i <= ndx;
|
||||
paradx = para;
|
||||
parady = para;
|
||||
paradx(seli) = paradx(seli) + dx(i)*dxscale(seli);
|
||||
parady(seli) = parady(seli) - dx(i)*dxscale(seli);
|
||||
paradxdy = paradx;
|
||||
paradxdy(seli) = paradxdy(seli) - dx(i)*dxscale(seli);
|
||||
fx = optmumlik20(para,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdx = optmumlik20(paradx,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdy = optmumlik20(parady,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdxdy = optmumlik20(paradxdy,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
gradx(i) = -( fx - fdx )/ (dx(i)*dxscale(seli));
|
||||
grady(i) = ( fx - fdy )/ (dx(i)*dxscale(seli));
|
||||
gradxy(i) = -(fx -fdxdy)/ sqrt( (dx(i)*dxscale(seli))^2 + (dx(i)*dxscale(seli))^2 );
|
||||
hessdiag(i) = -( 2*fx - fdx - fdy)/(dx(i)*dxscale(seli))^2;
|
||||
hessdiag(i) = -( fx - fdx - fdy + fdxdy )/(dx(i)*dx(i)*dxscale(seli)*dxscale(seli));
|
||||
i = i+1;
|
||||
end;
|
||||
% "Values";
|
||||
% -hessdiag';
|
||||
hessianfit(seli,seli) = -1*(hessdiag(1));
|
||||
% hessianfit(seli,seli) = -0.5*(hessdiag(3)+hessdiag(4));
|
||||
% locate 6,1;
|
||||
% "Value Used:" hessianfit[seli,seli];
|
||||
seli = seli+1;
|
||||
end;
|
||||
|
||||
%/* Now compute off-diagonal elements
|
||||
%** Make sure that correlations are between -1 and 1
|
||||
%** errorij contains the index of elements that are invalid
|
||||
%*/
|
||||
errorij = [ 0 0 0];
|
||||
|
||||
seli = 1;
|
||||
for seli = 1:npara;
|
||||
selj = seli+1
|
||||
while selj <= npara;
|
||||
% locate 1,1;
|
||||
% "hessianfit Element (" seli selj ")";
|
||||
i=1;
|
||||
while i <= ndx;
|
||||
paradx = para;
|
||||
parady = para;
|
||||
paradx(seli) = paradx(seli) + dx(i)*dxscale(seli);
|
||||
parady(selj) = parady(selj) - dx(i)*dxscale(selj);
|
||||
paradxdy = paradx;
|
||||
paradxdy(selj) = paradxdy(selj) - dx(i)*dxscale(selj);
|
||||
fx = optmumlik20(para,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdx = optmumlik20(paradx,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdy = optmumlik20(parady,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
fdxdy = optmumlik20(paradxdy,q,h2,aa,mf,gitno,gst,gobs,gend,gstart,nlags,rawdata,lpd);
|
||||
gradx(i) = -( fx - fdx )/ (dx(i)*dxscale(seli));
|
||||
grady(i) = ( fx - fdy )/ (dx(i)*dxscale(selj));
|
||||
gradxy(i) = -(fx -fdxdy)/ sqrt( (dx(i)*dxscale(selj))^2 + (dx(i)*dxscale(seli))^2 );
|
||||
hessdiag(i) = -( 2*fx - fdx - fdy)/(dx(i)*dxscale(seli))^2;
|
||||
hessdiag(i) = -( fx - fdx - fdy + fdxdy )/(dx(i)*dx(i)*dxscale(seli)*dxscale(selj));
|
||||
i = i+1
|
||||
end;
|
||||
% "Values";
|
||||
% -hessdiag';
|
||||
|
||||
% hessianfit(seli,selj) = -0.5*(hessdiag(3)+hessdiag(4));
|
||||
hessianfit(seli,selj) = -1*(hessdiag(1));
|
||||
|
||||
if ( hessianfit(seli,seli) == 0 ) or (hessianfit(selj,selj) == 0);
|
||||
corrij = 0;
|
||||
else;
|
||||
corrij = hessianfit(seli,selj)/sqrt(hessianfit(seli,seli)*hessianfit(selj,selj));
|
||||
end;
|
||||
|
||||
if (corrij < -1) or (corrij > 1);
|
||||
hessianfit(seli,selj)=0;
|
||||
errorij = [ errorij [seli selj corrij] ];
|
||||
end;
|
||||
hessianfit(selj,seli) = hessianfit(seli,selj);
|
||||
|
||||
% locate 6,1;
|
||||
% "Value Used: " hessianfit[seli,selj];
|
||||
% "Correlation:" corrij;
|
||||
% "Number of Errors:" rows(errorij)-1;
|
||||
selj=selj+1
|
||||
end;
|
||||
% seli = seli+1;
|
||||
end;
|
||||
|
||||
%cls;
|
||||
disp('Errors')
|
||||
disp(errorij);
|
||||
|
||||
|
||||
%/*******************************************************************************
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
%new;
|
||||
%closeall;
|
||||
%library user, pgraph, lbdlib;
|
||||
%format /mb1 /ros 16,8;
|
||||
%cls;
|
||||
|
||||
%/********************************************************
|
||||
%** Estimate the DSGE Model
|
||||
%** Models: 1 = RBC
|
||||
%** 2 = LBD
|
||||
%** 3 = LBD + Effort
|
||||
|
||||
%mspec = 3;
|
||||
%mprior = 2;
|
||||
%npara = 12;
|
||||
|
||||
%/********************************************************
|
||||
%** Import data on output growth and inflation: series (nobs,2)
|
||||
%** observations from 1954:III to 1997:IV:
|
||||
%**
|
||||
%** YY is composed of gdpq_cld and blsh_cl
|
||||
%*/
|
||||
%#include c:\projects\active\persist\Gauss\prog_t03\loaddata.g
|
||||
%loadm path = ^lpath para_names;
|
||||
%
|
||||
%/********************************************************
|
||||
%** Load Posterior Mode Estimates
|
||||
|
||||
%lpara = lpath $+ "\\" $+ lmodel $+ lprior $+ "mode";
|
||||
%open fhpara = ^lpara for read;
|
||||
%para = readr(fhpara,npara);
|
||||
%closeall fhpara;
|
||||
|
||||
%"Parameter | Estimate ";
|
||||
%outmat = para_names'~para;
|
||||
%let mask[1,2] = 0 1;
|
||||
%%let fmt[2,3] =
|
||||
% "-*.*s " 10 4
|
||||
% "*.*lf " 10 4;
|
||||
%d = printfm(outmat,(0 ~ 1),fmt);
|
||||
|
||||
%"";
|
||||
%"Prior*Likelihood at Mode";
|
||||
%fcn(para);
|
||||
%"Press Key to Continue";
|
||||
%k = keyw;
|
||||
%cls;
|
||||
% $$$
|
||||
% $$$ /*
|
||||
% $$$ goto evalhess;
|
||||
% $$$ */
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$ /* Initialize Output files
|
||||
% $$$ */
|
||||
% $$$ ohess = lpath $+ "\\" $+ lmodel $+ lprior $+ "hess";
|
||||
% $$$ create fhhess=^ohess with hess, npara, 8;
|
||||
% $$$ wr = writer(fhhess,hessianfit[1:npara,1:npara]);
|
||||
% $$$ closeall fhhess;
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$ /* Load hessianfit, compute penalty
|
||||
% $$$ */
|
||||
% $$$ evalhess:
|
||||
% $$$
|
||||
% $$$ lhess = lpath $+ "\\" $+ lmodel $+ lprior $+ "hess";
|
||||
% $$$ open fhhess = ^lhess for read;
|
||||
% $$$ HHm = readr( fhhess,npara );
|
||||
% $$$ closeall fhhess;
|
||||
% $$$
|
||||
% $$$ /* hessianfit is of reduced rank. Keep zero rows and columns and do SVD
|
||||
% $$$ */
|
||||
% $$$ if mspec == 1;
|
||||
% $$$ rankHHm = 9;
|
||||
% $$$ elseif mspec == 2;
|
||||
% $$$ rankHHm = 11;
|
||||
% $$$ else;
|
||||
% $$$ rankHHm = 12;
|
||||
% $$$ endif;
|
||||
% $$$
|
||||
% $$$ /* Create Inverse by Singular Value Decomposition
|
||||
% $$$ */
|
||||
% $$$ {u , s, v} = svd1(HHM);
|
||||
% $$$ invHHMdet = 1;
|
||||
% $$$
|
||||
% $$$ i = 1;
|
||||
% $$$ do until i > npara;
|
||||
% $$$ if i > rankHHM;
|
||||
% $$$ s[i,i] = 0;
|
||||
% $$$ else;
|
||||
% $$$ s[i,i] = 1/s[i,i];
|
||||
% $$$ invHHMdet = invHHMdet*s[i,i];
|
||||
% $$$ endif;
|
||||
% $$$ i = i+1;
|
||||
% $$$ endo;
|
||||
% $$$
|
||||
% $$$ invHHM = u*s*u';
|
||||
% $$$ sigmult = u*sqrt(s);
|
||||
% $$$
|
||||
% $$$ "Determinant of minus hessianfit";
|
||||
% $$$ invHHMdet;
|
||||
% $$$ "sqrt(Diagonal of Inverse hessianfit)";
|
||||
% $$$ sqrt(diag(invHHM));
|
||||
% $$$
|
||||
% $$$ "Post Mode Penalty";
|
||||
% $$$ penalt = (rankHHM/2)*ln(2*pi) + 0.5*ln(invHHMdet);
|
||||
% $$$ penalt;
|
||||
% $$$
|
||||
% $$$ /* Initialize Output files
|
||||
% $$$ */
|
||||
% $$$ omult = lpath $+ "\\" $+ lmodel $+ lprior $+ "mult";
|
||||
% $$$ create fhmult=^omult with MULT, npara, 8;
|
||||
% $$$ wr = writer(fhmult,sigmult);
|
||||
% $$$
|
||||
% $$$ closeall fhmult;
|
||||
% $$$ end;
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$ %/****************************************************/
|
||||
% $$$ %/* PROCEDURES */
|
||||
% $$$ %/****************************************************/
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$ %proc (1) = fcn(para);
|
||||
% $$$ %local lnpY, lnprio1, lnprio2, obsmean, obsvar;
|
||||
% $$$
|
||||
% $$$ %{lnpy, obsmean, obsvar} = evallbd( para,mspec,T0,YY);
|
||||
% $$$
|
||||
% $$$ %/* Evaluate the Prior density
|
||||
% $$$
|
||||
% $$$ % lnprio1 = priodens( para, pmean, pstdd, pshape);
|
||||
% $$$ % lnprio2 = priomuphi( para );
|
||||
% $$$
|
||||
% $$$ %retp(real(lnpY+lnprio1+lnprio2));
|
||||
% $$$ %endp;
|
||||
% $$$
|
||||
% $$$ /***************************************************************************
|
||||
% $$$ */
|
||||
% $$$
|
||||
% $$$ proc (1) = priomuphi(para);
|
||||
% $$$ local muphi, lnprio;
|
||||
% $$$ muphi = para[7:8];
|
||||
% $$$ if mspec > 1;
|
||||
% $$$ lnprio = -ln(2*pi) - 0.5*ln(det(muphi_v0))
|
||||
% $$$ - 0.5*(muphi - muphi_m0)'*inv(muphi_v0)*(muphi - muphi_m0);
|
||||
% $$$ else;
|
||||
% $$$ lnprio = 0;
|
||||
% $$$ endif;
|
||||
% $$$ retp(lnprio);
|
||||
% $$$ endp;
|
||||
% $$$
|
||||
% $$$
|
||||
% $$$
|
||||
|
|
@ -8,8 +8,8 @@ function make_ex_
|
|||
if isempty(oo_.exo_steady_state)
|
||||
oo_.exo_steady_state = zeros(M_.exo_nbr,1);
|
||||
end
|
||||
if M_.exo_det_nbr > 1 & isempty(oo_.exo_det_steadystate)
|
||||
oo_.exo_det_steadystate = zeros(M_.exo_det_nbr,1);
|
||||
if M_.exo_det_nbr > 1 & isempty(oo_.exo_det_steady_state)
|
||||
oo_.exo_det_steady_state = zeros(M_.exo_det_nbr,1);
|
||||
end
|
||||
if isempty(oo_.exo_simul)
|
||||
if isempty(ex0_)
|
||||
|
|
@ -36,24 +36,24 @@ function make_ex_
|
|||
if M_.exo_det_nbr > 0
|
||||
if isempty(oo_.exo_det_simul)
|
||||
if isempty(ex_det0_)
|
||||
oo_.exo_det_simul = [ones(ykmin_+options_.periods+ykmax_,1)*M_.exo_det_steadystate'];
|
||||
oo_.exo_det_simul = [ones(M_.maximum_lag+options_.periods+M_.maximum_lead,1)*oo_.exo_det_steady_state'];
|
||||
else
|
||||
oo_.exo_det_simul = [ones(ykmin_,1)*ex_det0_';ones(options_.periods+ykmax_,1)*M_.exo_det_steadystate'];
|
||||
oo_.exo_det_simul = [ones(M_.maximum_lag,1)*ex_det0_';ones(options_.periods+M_.maximum_lead,1)*oo_.exo_det_steady_state'];
|
||||
end
|
||||
elseif size(oo_.exo_det_simul,2) < length(M_.exo_det_steadystate)
|
||||
k = size(oo_.exo_det_simul,2)+1:length(M_.exo_det_steadystate);
|
||||
elseif size(oo_.exo_det_simul,2) < length(oo_.exo_det_steady_state)
|
||||
k = size(oo_.exo_det_simul,2)+1:length(oo_.exo_det_steady_state);
|
||||
if isempty(ex_det0_)
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul ones(ykmin_+size(oo_.exo_det_simul,1)+ykmax_,1)*M_.exo_det_steadystate(k)'];
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul ones(M_.maximum_lag+size(oo_.exo_det_simul,1)+M_.maximum_lead,1)*oo_.exo_det_steady_state(k)'];
|
||||
else
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul [ones(ykmin_,1)*ex_det0_(k)'; ones(size(oo_.exo_det_simul,1)-ykmin_+ykmax_, ...
|
||||
1)*M_.exo_det_steadystate(k)']];
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul [ones(M_.maximum_lag,1)*ex_det0_(k)'; ones(size(oo_.exo_det_simul,1)-M_.maximum_lag+M_.maximum_lead, ...
|
||||
1)*oo_.exo_det_steady_state(k)']];
|
||||
end
|
||||
elseif size(oo_.exo_det_simul,1) < ykmin_+ykmax_+options_.periods
|
||||
elseif size(oo_.exo_det_simul,1) < M_.maximum_lag+M_.maximum_lead+options_.periods
|
||||
if isempty(ex_det0_)
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul; ones(ykmin_+options_.periods+ykmax_-size(oo_.exo_det_simul,1),1)*M_.exo_det_steadystate'];
|
||||
oo_.exo_det_simul = [oo_.exo_det_simul; ones(M_.maximum_lag+options_.periods+M_.maximum_lead-size(oo_.exo_det_simul,1),1)*oo_.exo_det_steady_state'];
|
||||
else
|
||||
oo_.exo_det_simul = [ones(ykmin_,1)*ex_det0_'; oo_.exo_det_simul; ones(options_.periods+ykmax_-size(oo_.exo_det_simul, ...
|
||||
1),1)*M_.exo_det_steadystate'];
|
||||
oo_.exo_det_simul = [ones(M_.maximum_lag,1)*ex_det0_'; oo_.exo_det_simul; ones(options_.periods+M_.maximum_lead-size(oo_.exo_det_simul, ...
|
||||
1),1)*oo_.exo_det_steady_state'];
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,226 +0,0 @@
|
|||
% Program calculating the posterior density
|
||||
% 1. define xparam
|
||||
% 2. call model setup & reduction program
|
||||
% 3. prepare state space variables and kalman-filter setup
|
||||
% 4. evaluate likelihood with kalman filter
|
||||
% 5. evaluate prior
|
||||
%------------------------------------------------------------------------------
|
||||
%------------------------------------------------------------------------------
|
||||
function [fval,cost_flag,atT,innov,ys,trend_coeff] = ...
|
||||
mj_optmumlik(xparam1,gend,rawdata,algo);
|
||||
|
||||
% algo = 1: computes filter + likelihood
|
||||
% alog = 2: computes filter + likelihood + smoother
|
||||
|
||||
global M_ oo_ estim_params_ options_
|
||||
global bayestopt_ dr1_test_ trend_coeff_ xparam_test
|
||||
|
||||
xparam_test = xparam1;
|
||||
cost_flag = 1;
|
||||
if options_.mode_compute ~= 1 && any(xparam1 < bayestopt_.lb)
|
||||
k = find(xparam1 < bayestopt_.lb);
|
||||
fval = bayestopt_.penalty+sum(bayestopt_.lb(k)-xparam1(k));
|
||||
cost_flag = 0;
|
||||
return;
|
||||
end
|
||||
if options_.mode_compute ~= 1 && any(xparam1 > bayestopt_.ub)
|
||||
k = find(xparam1 > bayestopt_.ub);
|
||||
fval = bayestopt_.penalty+sum(xparam1(k)-bayestopt_.ub(k));
|
||||
cost_flag = 0;
|
||||
return;
|
||||
end
|
||||
|
||||
nobs = size(options_.varobs,1);
|
||||
|
||||
q = M_.Sigma_e;
|
||||
for i=1:estim_params_.nvx
|
||||
k =estim_params_.var_exo(i,1);
|
||||
q(k,k) = xparam1(i)*xparam1(i);
|
||||
end
|
||||
|
||||
offset = estim_params_.nvx;
|
||||
h = zeros(nobs,nobs);
|
||||
for i=1:estim_params_.nvn
|
||||
k =estim_params_.var_endo(i,1);
|
||||
h(k,k) = xparam1(i+offset)*xparam1(i+offset);
|
||||
end
|
||||
|
||||
offset = offset+estim_params_.nvn;
|
||||
for i=1:estim_params_.ncx
|
||||
k1 =estim_params_.corrx(i,1);
|
||||
k2 =estim_params_.corrx(i,2);
|
||||
q(k1,k2) = xparam1(i+offset)*sqrt(q(k1,k1)*q(k2,k2));
|
||||
q(k2,k1) = q(k1,k2);
|
||||
end
|
||||
|
||||
offset = offset+estim_params_.ncx;
|
||||
for i=1:estim_params_.ncn
|
||||
k1 =estim_params_.corrn(i,1);
|
||||
k2 =estim_params_.corrn(i,2);
|
||||
h(k1,k2) = xparam1(i+offset)*sqrt(h(k1,k1)*h(k2,k2));
|
||||
h(k2,k1) = h(k1,k2);
|
||||
end
|
||||
|
||||
offset = offset+estim_params_.ncn;
|
||||
for i=1:estim_params_.np
|
||||
%this dosn't work with struct field M_.param
|
||||
%assignin('base',deblank(estim_params_.param_names(i,:)),xparam1(i+offset));
|
||||
tmpvalue = xparam1(i+offset);
|
||||
eval('caller',[ deblank(estim_params_.param_names(i,:)) ' = tmpvalue;']);
|
||||
end
|
||||
|
||||
%------------------------------------------------------------------------------
|
||||
% 2. call model setup & reduction program
|
||||
%------------------------------------------------------------------------------
|
||||
|
||||
[A,B,ys] = dynare_resolve;
|
||||
|
||||
if dr1_test_(1) == 1
|
||||
fval = bayestopt_.penalty*exp(dr1_test_(2));
|
||||
cost_flag = 0;
|
||||
return;
|
||||
elseif dr1_test_(1) == 2;
|
||||
fval = bayestopt_.penalty*exp(dr1_test_(2));
|
||||
cost_flag = 0;
|
||||
return;
|
||||
elseif dr1_test_(1) == 3;
|
||||
fval = bayestopt_.penalty*exp(dr1_test_(2));
|
||||
cost_flag = 0;
|
||||
return;
|
||||
end
|
||||
|
||||
if options_.loglinear == 1
|
||||
ys1 = log(ys(bayestopt_.mfys));
|
||||
else
|
||||
ys1 = ys(bayestopt_.mfys);
|
||||
end
|
||||
|
||||
aa = A;
|
||||
|
||||
np = size(A,1);
|
||||
mf = eye(np);
|
||||
mf = bayestopt_.mf;
|
||||
|
||||
% Set initial values @
|
||||
|
||||
|
||||
at = zeros(np,gend+1);
|
||||
gconst = log(2*pi);
|
||||
lik = zeros(gend,1);
|
||||
|
||||
if options_.lik_init == 1
|
||||
p0 = lyapunov_symm(aa,B*q*B');
|
||||
elseif options_.lik_init == 2
|
||||
p0=eye(np)*10.0;
|
||||
end
|
||||
pt = p0;
|
||||
BqB = B*q*B';
|
||||
|
||||
trend_coeff = zeros(nobs,1);
|
||||
if bayestopt_.with_trend == 1
|
||||
nx1 = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+ ...
|
||||
estim_params_.ncn;
|
||||
for i=1:nobs
|
||||
trend_coeff(i) = eval(bayestopt_.trend_coeff{i});
|
||||
end
|
||||
end
|
||||
|
||||
not_steady = 1;
|
||||
ldetf_old = NaN;
|
||||
warning_state = warning;
|
||||
warning off;
|
||||
if algo == 1
|
||||
for t = 1:gend
|
||||
if not_steady
|
||||
ptt1 = aa*pt*aa'+BqB;
|
||||
f = ptt1(mf,mf)+h;
|
||||
ldetf = log(det(f));
|
||||
finv = inv(f);
|
||||
if any(isinf(finv)) | ~isreal(ldetf)
|
||||
disp('singularity in Kalman filter');
|
||||
fval = bayestopt_.penalty;
|
||||
warning(warning_state);
|
||||
cost_flag = 0;
|
||||
dr1_test_(1) = 4;
|
||||
return
|
||||
end
|
||||
pt = ptt1-ptt1(:,mf)*finv*ptt1(mf,:);
|
||||
if abs(ldetf-ldetf_old) < 1e-12
|
||||
not_steady = 0;
|
||||
end
|
||||
ldetf_old = ldetf;
|
||||
end
|
||||
att1 = aa*at(:,t);
|
||||
v = rawdata(t,:)'-att1(mf,:)-ys1;
|
||||
if bayestopt_.with_trend == 1
|
||||
v = v - trend_coeff*t;
|
||||
end
|
||||
at(:,t+1) = att1+ptt1(:,mf)*finv*v;
|
||||
if t > options_.presample
|
||||
lik(t,1) = ldetf+v'*finv*v;
|
||||
end
|
||||
end
|
||||
elseif algo == 2
|
||||
PPt = zeros(np^2,gend);
|
||||
for t = 1:gend
|
||||
if not_steady
|
||||
ptt1 = aa*pt*aa'+BqB;
|
||||
f = ptt1(mf,mf)+h;
|
||||
ldetf = log(det(f));
|
||||
finv = inv(f);
|
||||
if any(isinf(finv)) | ~isreal(ldetf)
|
||||
% disp('singularity in Kalman filter');
|
||||
fval = bayestopt_.penalty;
|
||||
warning(warning_state);
|
||||
cost_flag = 0;
|
||||
dr1_test_(1) = 4;
|
||||
return
|
||||
end
|
||||
pt = ptt1-ptt1(:,mf)*finv*ptt1(mf,:);
|
||||
if abs(ldetf-ldetf_old) < 1e-12
|
||||
not_steady = 0;
|
||||
end
|
||||
ldetf_old = ldetf;
|
||||
end
|
||||
att1 = aa*at(:,t);
|
||||
v = rawdata(t,:)' - att1(mf,:) - ys1;
|
||||
if bayestopt_.with_trend == 1
|
||||
v = v - trend_coeff*t;
|
||||
end
|
||||
at(:,t+1) = att1+ptt1(:,mf)*finv*v;
|
||||
PPt(:,t) = pt(:);
|
||||
if t > options_.presample
|
||||
lik(t,1) = ldetf+v'*finv*v;
|
||||
end
|
||||
end
|
||||
atT =zeros(np,gend+1);
|
||||
atT(:,gend+1) = at(:,gend+1);
|
||||
innov = zeros(M_.exo_nbr,gend);
|
||||
for t = gend:-1:1
|
||||
pt = reshape(PPt(:,t),np,np);
|
||||
ptt1 = aa*pt*aa'+BqB;
|
||||
Pstar = pt*aa'*pinv(ptt1);
|
||||
atT(:,t) = at(:,t)+Pstar*(atT(:,t+1)-aa*at(:,t));
|
||||
shocks1 = atT(:,t+1)-aa*atT(:,t);
|
||||
innov(:,t) = B\shocks1;
|
||||
end
|
||||
end
|
||||
|
||||
warning(warning_state);
|
||||
|
||||
likelihood = 0.5*sum(nobs*gconst+lik(options_.presample+1:end));
|
||||
|
||||
if imag(likelihood) ~= 0
|
||||
|
||||
likelihood = 10000000;
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------------------
|
||||
% PRIOR SPECIFICATION
|
||||
% ------------------------------------------------------------------------------
|
||||
|
||||
lnprior = priordens(xparam1, bayestopt_.pshape, bayestopt_.p1, bayestopt_.p2, bayestopt_.p3, bayestopt_.p4 );
|
||||
|
||||
fval = (likelihood-lnprior);
|
||||
|
||||
% 11/18/03 MJ changed input parameters for priordens()
|
||||
|
|
@ -11,6 +11,7 @@ global it_
|
|||
|
||||
|
||||
options_ = set_default_option(options_,'olr',0);
|
||||
options_ = set_default_option(options_,'jacobian_flag',1);
|
||||
info = 0;
|
||||
|
||||
it_ = M_.maximum_lag + 1 ;
|
||||
|
|
@ -24,17 +25,18 @@ tempex = oo_.exo_simul;
|
|||
oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+M_.maximum_lead+1,1);
|
||||
if M_.exo_det_nbr > 0
|
||||
tempexdet = oo_.exo_det_simul;
|
||||
oo_.exo_det_simul = ones(M_.maximum_lag+1,1)*oo_.exo_steady_statedet_';
|
||||
oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+M_.maximum_lead+1,1);
|
||||
end
|
||||
dr.ys = ys;
|
||||
fh = str2func([M_.fname '_static']);
|
||||
if options_.linear == 0
|
||||
if max(abs(feval(fh,dr.ys,oo_.exo_steady_state))) > options_.dynatol & options_.olr == 0
|
||||
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state]))) > options_.dynatol & options_.olr == 0
|
||||
if exist([M_.fname '_steadystate'])
|
||||
[dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,oo_.exo_steady_state);
|
||||
[dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,...
|
||||
[oo_.exo_steady_state; oo_.exo_det_steady_state]);
|
||||
else
|
||||
%[dr.ys,check1] = dynare_solve(fh,dr.ys,oo_.exo_steady_state);
|
||||
[dr.ys,check1] = dynare_solve(fh,dr.ys,1,oo_.exo_steady_state);
|
||||
[dr.ys,check1] = dynare_solve(fh,dr.ys,options_.jacobian_flag,...
|
||||
[oo_.exo_steady_state; oo_.exo_det_steady_state]);
|
||||
end
|
||||
if check1
|
||||
info(1) = 20;
|
||||
|
|
|
|||
|
|
@ -4,6 +4,8 @@ function steady(linear)
|
|||
|
||||
global M_ oo_ options_ ys0_
|
||||
|
||||
options_ = set_default_option(options_,'jacobian_flag',1);
|
||||
|
||||
steady_;
|
||||
|
||||
disp(' ')
|
||||
|
|
|
|||
|
|
@ -2,23 +2,19 @@
|
|||
%
|
||||
function steady_()
|
||||
|
||||
global M_ oo_ it_
|
||||
|
||||
|
||||
x = oo_.steady_state ;
|
||||
xlen = M_.maximum_lag + M_.maximum_lead + 1 ;
|
||||
nn = size(M_.lead_lag_incidence,2) ;
|
||||
it_ = M_.maximum_lag+1 ;
|
||||
x = repmat(oo_.exo_steady_state',xlen,1);
|
||||
|
||||
if M_.exo_det_nbr > 0
|
||||
x = [x, repmat(oo_.exo_det_steady_state',M_.maximum_lag+1,1)] ;
|
||||
end
|
||||
global M_ oo_ it_ options_
|
||||
|
||||
if exist([M_.fname '_steadystate'])
|
||||
[oo_.steady_state,check] = feval([M_.fname '_steadystate'],oo_.steady_state,x);
|
||||
[oo_.steady_state,check] = feval([M_.fname '_steadystate'],...
|
||||
oo_.steady_state,...
|
||||
[oo_.exo_steady_state; ...
|
||||
oo_.exo_det_steady_state]);
|
||||
else
|
||||
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],oo_.steady_state,1,x);
|
||||
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
|
||||
oo_.steady_state,...
|
||||
options_.jacobian_flag, ...
|
||||
[oo_.exo_steady_state; ...
|
||||
oo_.exo_det_steady_state]);
|
||||
end
|
||||
|
||||
if check ~= 0
|
||||
|
|
|
|||
Loading…
Reference in New Issue