diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m deleted file mode 100644 index 6cb1c3101..000000000 --- a/matlab/DsgeLikelihood_hh.m +++ /dev/null @@ -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)]; diff --git a/matlab/check_model.m b/matlab/check_model.m index 7189ca248..c7e261043 100644 --- a/matlab/check_model.m +++ b/matlab/check_model.m @@ -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 diff --git a/matlab/dr1.m b/matlab/dr1.m index d6c159e16..0a2fb2b1b 100644 --- a/matlab/dr1.m +++ b/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); diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe index 3fbf7fe8a..eaada2ce5 100755 Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ diff --git a/matlab/hessian2.m b/matlab/hessian2.m deleted file mode 100644 index e113b8f8c..000000000 --- a/matlab/hessian2.m +++ /dev/null @@ -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; -% $$$ -% $$$ -% $$$ diff --git a/matlab/make_ex_.m b/matlab/make_ex_.m index 2122d24e0..c20878acd 100644 --- a/matlab/make_ex_.m +++ b/matlab/make_ex_.m @@ -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 diff --git a/matlab/mj_optmumlik.m b/matlab/mj_optmumlik.m deleted file mode 100644 index f7f4f8dfb..000000000 --- a/matlab/mj_optmumlik.m +++ /dev/null @@ -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() \ No newline at end of file diff --git a/matlab/resol.m b/matlab/resol.m index c54df7bf0..10dc3e0b5 100644 --- a/matlab/resol.m +++ b/matlab/resol.m @@ -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; diff --git a/matlab/steady.m b/matlab/steady.m index e18b79896..5bb063f9c 100644 --- a/matlab/steady.m +++ b/matlab/steady.m @@ -4,6 +4,8 @@ function steady(linear) global M_ oo_ options_ ys0_ + options_ = set_default_option(options_,'jacobian_flag',1); + steady_; disp(' ') diff --git a/matlab/steady_.m b/matlab/steady_.m index 0ad3e296f..3b5326226 100644 --- a/matlab/steady_.m +++ b/matlab/steady_.m @@ -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