diff --git a/matlab/CreateBenchmark.m b/matlab/CreateBenchmark.m deleted file mode 100644 index 592810ab0..000000000 --- a/matlab/CreateBenchmark.m +++ /dev/null @@ -1,8 +0,0 @@ -function CreateBenchmark - -% stephane.adjemian@cepremap.cnrs.fr [12-06-2004] - -global oo_ - -eval([M_.fname '_oo_ = oo_;']) -eval(['save ' M_.fname '_benchmark_oo ' M_.fname '_oo_']) \ No newline at end of file diff --git a/matlab/InitSeed.m b/matlab/InitSeed.m deleted file mode 100644 index cf3c0824c..000000000 --- a/matlab/InitSeed.m +++ /dev/null @@ -1,45 +0,0 @@ -function InitSeed - -% stephane.adjemian@cepremap.cnrs.fr [12-02-2004] - -normalseed = [ 220271 ; - 220378]; - -uniformseed = [ 0.6776 ; ... - 0.8239 ; ... - 0.6219 ; ... - 0.8169 ; ... - 0.2796 ; ... - 0.4337 ; ... - 0.6205 ; ... - 0.8994 ; ... - 0.9718 ; ... - 0.7993 ; ... - 0.4444 ; ... - 0.8896 ; ... - 0.9231 ; ... - 0.5838 ; ... - 0.1788 ; ... - 0.6556 ; ... - 0.7464 ; ... - 0.7545 ; ... - 0.2400 ; ... - 0.4431 ; ... - 0.3453 ; ... - 0.6687 ; ... - 0.6287 ; ... - 0.8371 ; ... - 0.8041 ; ... - 0.6802 ; ... - 0.5864 ; ... - 0.7713 ; ... - 0.8282 ; ... - 0.5748 ; ... - 0.0999 ; ... - 0.6360 ; ... - 0 ; ... - 0.0000 ; ... - 0.0000]; - -randn('state',normalseed); -rand('state',uniformseed); \ No newline at end of file diff --git a/matlab/check_mh.m b/matlab/check_mh.m deleted file mode 100644 index 883641a84..000000000 --- a/matlab/check_mh.m +++ /dev/null @@ -1,24 +0,0 @@ -function check_mh(fname) - eval(['load ' fname]); - nb = size(x3,3); - np = size(x2,2); - nr = size(x2,1); - - j1 = ceil(0.5*nr); - x = [j1:100:nr]; - z = []; - for i=1:np - y1 = zeros(size(x),nb); - for k=1:nb - for j=1:length(x) - y1(j,k) = mean(x2(1:x(j),i,k)); - end - end - my_subplot(i,np,4,5,'MH convergence'); - plot([y1]) - xmin = min(min(x2(:,i,:))); - xmax = max(max(x2(:,i,:))); - z = [z; [i sum(sum(x3(:,i,:) < xmin))/nr sum(sum(x3(:,i) > xmax))/nr]]; - end - disp(z) - \ No newline at end of file diff --git a/matlab/compDist.m b/matlab/compDist.m deleted file mode 100644 index 555015188..000000000 --- a/matlab/compDist.m +++ /dev/null @@ -1,171 +0,0 @@ -function compdist(xparam1, x2, pltopt, figurename) -global bayestopt_ estim_params_ M_ options_ - -% NOTE: If pltopt ~= 'All' compdist.m just draws prior densities. - -%% Set density estimation parameters: -number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two. -bandwidth = 0; % Rule of thumb optimal bandwidth parameter. -kernel_function = 'gaussian'; % You can switch to: 'epanechnikov', 'quartic', 'triangle', -% 'triweight', 'uniform' or 'cosinus' kernels (iff bandwidth=0, see posterior_density_estimate.m). -truncprior = 10^(-3); - - -npar=length(xparam1); -nruns=length(x2); -icol=ceil(sqrt(npar)); -iraw=icol; -if (icol-1)*(icol-2)>=npar - iraw = icol-2; - icol=icol-1; -elseif (icol)*(icol-2)>=npar - iraw = icol-2; -elseif icol*(icol-1)>=npar - iraw=icol-1; -end - -pmean=bayestopt_.pmean; -pshape=bayestopt_.pshape; -p1 = bayestopt_.p1; -p2 = bayestopt_.p2; -p3 = bayestopt_.p3; -p4 = bayestopt_.p4; - -figure('Name',figurename) -for i=1:npar; - if i<=estim_params_.nvx - vname = deblank(M_.exo_name(estim_params_.var_exo(i,1),:)); - nam=['SE_{',vname,'}']; - elseif i<=(estim_params_.nvx+estim_params_.nvn) - deblank(options_.varobs(estim_params_.var_endo(i-estim_params_.nvx,1),:)); - nam=['SE_{EOBS_',vname,'}']; - elseif i<=(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx) - j = i - (estim_params_.nvx+estim_params_.nvn); - k1 = estim_params_.corrx(j,1); - k2 = estim_params_.corrx(j,2); - vname = [deblank(M_.exo_name(k1,:)) ',' deblank(M_.exo_name(k2,:))]; - nam=['CC_{',vname,'}']; - elseif i<=(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+ ... - estim_params_.ncn) - j = i - (estim_params_.nvx+estim_params_.nvn+estim_params_.ncx); - k1 = estim_params_.corrn(j,1); - k2 = estim_params_.corrn(j,2); - vname = [deblank(M_.exo_name(k1,:)) ',' deblank(M_.exo_name(k2,:))]; - nam=['CC_{EOBS_',vname,'}']; - else - j = i - (estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+ ... - estim_params_.ncn); - nam = deblank(estim_params_.param_names(j,:)); - end - subplot(iraw, icol, i); - if strcmpi(pltopt,'all'); % Estimation of the density... - [abscissa,ff,h] = posterior_density_estimate(x2(round(options_.mh_drop*nruns):end,i),... - number_of_grid_points,bandwidth,kernel_function); - plot(abscissa,ff,'-k','linewidth',2); - end; - a = 0; - b = 0; - if pshape(i) == 1; %/* BETA Prior */ - density = inline('((1-x).^(b-1)).*x.^(a-1)./beta(a,b)','x','a','b'); - mu = (p1(i)-p3(i))/(p4(i)-p3(i)); - stdd = p2(i)/(p4(i)-p3(i)); - a = (1-mu)*mu^2/stdd^2 - mu; - b = a*(1/mu - 1); - infbound = qbeta(truncprior,a,b); - supbound = qbeta(1-truncprior,a,b); - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = density(abscissa,a,b); - abscissa = abscissa*(p4(i)-p3(i))+p3(i); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - elseif pshape(i) == 2; %/* GAMMA PRIOR */ -% density = -% inline('((x/b).^(a-1)).*exp(-x/b)*inv(b*gamma(a))','x','a','b'); - mu = p1(i)-p3(i); - b = p2(i)^2/mu; - a = mu/b; - infbound = mj_qgamma(truncprior,a)*b; - supbound = mj_qgamma(1-truncprior,a)*b; - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = exp(lpdfgam(abscissa,a,b)); - abscissa = abscissa + p3(i); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - elseif pshape(i) == 3; %/* GAUSSIAN PRIOR */ - density = inline('inv(sqrt(2*pi)*b)*exp(-0.5*((x-a)/b).^2)','x','a','b'); - a = p1(i); - b = p2(i); - infbound = qnorm(truncprior,a,b); - supbound = qnorm(1-truncprior,a,b); - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = density(abscissa,a,b); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - elseif pshape(i) == 4; %/* INVGAMMA PRIOR type 1 */ - density = inline('2*inv(gamma(nu/2))*(x.^(-nu-1))*((s/2)^(nu/2)).*exp(-s./(2*x.^2))','x','s','nu'); - nu = p2(i); - s = p1(i); - a = nu/2; - b = 2/s; - infbound = 1/sqrt(mj_qgamma(1-10*truncprior,a)*b); - supbound = 1/sqrt(mj_qgamma(10*truncprior,a)*b); - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = density(abscissa,s,nu); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - elseif pshape(i) == 5; %/* UNIFORM PRIOR */ - density = inline('(x.^0)/(b-a)','x','a','b'); - a = p1(i); - b = p2(i); - infbound = a; - supbound = b; - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = density(abscissa,a,b); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - elseif pshape(i) == 6; %/* INVGAMMA PRIOR type 2 */ - density = inline('inv(gamma(nu/2))*(x.^(-.5(nu+2)))*((s/2)^(nu/2)).*exp(-s./(2*x))','x','s','nu'); - nu = p2(i); - s = p1(i); - a = nu/2; - b = 2/s; - infbound = 1/(qgamma(1-truncprior,a)*b); - supbound = 1/(qgamma(truncprior,a)*b); - stepsize = (supbound-infbound)/200; - abscissa = infbound:stepsize:supbound; - f = density(abscissa,s,nu); - if strcmp(pltopt,'all'); - top = max([max(ff);max(f)]); - end; - end; - hold on; - k = [1:length(f)]; - if pshape(i) ~= 5 - [junk,k1] = max(f); - if k1 == 1 | k1 == length(f) - k = find(f < 10); - end - end - hh = plot(abscissa(k),f(k),'-k','linewidth',2); - set(hh,'color',[0.7 0.7 0.7]); - %hold on; %% ?! - if strcmp(pltopt,'all'); - plot( [xparam1(i) xparam1(i)], [0,top], '--g', 'linewidth', 2); - end; - title(nam,'Interpreter','none'); - hold off; -end; -drawnow - -% 12/01/03 MJ adapted from M. Ratto's version - diff --git a/matlab/dlognorm.m b/matlab/dlognorm.m deleted file mode 100644 index 4dcfed550..000000000 --- a/matlab/dlognorm.m +++ /dev/null @@ -1,16 +0,0 @@ -function f = dlognorm(x,lambda,zeta) -%DLOGNORM The log-normal density function -% -% f = dlognorm(x,lambda,zeta) - -% Copyright (c) Halfdan Grage, Anders Holtsberg - -if any(any(zeta<=0)) - error('Parameter zeta is wrong') -end - -neg = find(x <= 0); -[n1,n2] = size(neg); -x(neg) = ones(n1,n2); -f = dnorm(log(x),lambda,zeta)./x; -f(neg) = zeros(n1,n2); diff --git a/matlab/dr2.m b/matlab/dr2.m deleted file mode 100644 index 675513cdc..000000000 --- a/matlab/dr2.m +++ /dev/null @@ -1,11 +0,0 @@ -% Copyright (C) 2001 Michel Juillard -% -% function used for dr_algo == 1 -function ghs2=dr2(ys,dr) - global M_ - - dr.ys = ys; - fh = str2func([M_.fname '_static']); - dr.fbias = 2*feval(fh,dr.ys); - dr=dr1(dr,0); - ghs2 = dr.ghs2; diff --git a/matlab/fnorm.m b/matlab/fnorm.m deleted file mode 100644 index 2c99b2b15..000000000 --- a/matlab/fnorm.m +++ /dev/null @@ -1,2 +0,0 @@ -function e = fnorm(p2,p,p1,perc) - e = p - pnorm(perc,p1,p2); \ No newline at end of file diff --git a/matlab/lnsrch.m b/matlab/lnsrch.m deleted file mode 100644 index 80e9cdae3..000000000 --- a/matlab/lnsrch.m +++ /dev/null @@ -1,86 +0,0 @@ -% Copyright (C) 2001 Michel Juillard -% -function [x,f,fvec,check]=lnsrch(xold,fold,g,p,stpmax,func,varargin) - - alf = 1e-4 ; - tolx = 3.7e-11 ; - alam = 1; - - nn = size(xold,1) ; - summ = sqrt(sum(p.*p)) ; - if ~isfinite(summ) - error(['Some element of Newton direction isn''t finite. Jacobian maybe' ... - ' singular or there is a problem with initial values']) - end - - if summ > stpmax - p=p.*stpmax/summ ; - end - - slope = g'*p ; - - test = max(abs(p)'./max([abs(xold)';ones(1,nn)])) ; - alamin = tolx/test ; - - if alamin > 0.1 - alamin = 0.1; - end - - while 1 - if alam < alamin - check = 1 ; - return - end - - x = xold + (alam*p) ; - fvec = feval(func,x,varargin{:}) ; - f = 0.5*fvec'*fvec ; - - if any(isnan(fvec)) - alam = alam/2 ; - alam2 = alam ; - f2 = f ; - fold2 = fold ; - else - - if f <= fold+alf*alam*slope - check = 0; - break ; - else - if alam == 1 - 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 - tmplam = -slope/(2*b) ; - else - disc = (b^2)-3*a*slope ; - - if disc < 0 - error ('Roundoff problem in nlsearch') ; - else - tmplam = (-b+sqrt(disc))/(3*a) ; - end - - end - - if tmplam > 0.5*alam - tmplam = 0.5*alam; - end - - end - - alam2 = alam ; - f2 = f ; - fold2 = fold ; - alam = max([tmplam;(0.1*alam)]) ; - end - end - end - -% 01/14/01 MJ lnsearch is now a separate function -% 01/12/03 MJ check for finite summ to avoid infinite loop when Jacobian -% is singular or model is denormalized \ No newline at end of file diff --git a/matlab/p2toperc.m b/matlab/p2toperc.m deleted file mode 100644 index fc6cb0c74..000000000 --- a/matlab/p2toperc.m +++ /dev/null @@ -1,29 +0,0 @@ -function y=p2toperc(x) - - n = size(x,1); - y = zeros(n,1); - - for i=1:n -% if x(i,6) == 0 have to wait until the 2nd part works - if 1 - y(i) = x(i,3); - else - if x(i,3) < x(i,2) - p = 0.05; - else - p = 0.95; - end - - if x(i,1) == 1 - y(i) = fsolve(@fbeta,1,p,x(i,2),x(i,3)); - elseif x(i,1) == 2 - y(i) = fsolve(@fgamma,1,p,x(i,2),x(i,3)); - elseif x(i,1) == 3 - y(i) = (x(i,3)-x(i,2))/qnorm(p,0,1); - elseif x(i,1) == 4 - y(i) = fsolve(@figamm,1,p,x(i,2),x(i,3)); - elseif x(i,1) == 5 -% y(i) = fsolve(@fgamma,1,p,x(i,2)); - end - end - end \ No newline at end of file diff --git a/matlab/resol1.m b/matlab/resol1.m deleted file mode 100644 index 1646fcbcc..000000000 --- a/matlab/resol1.m +++ /dev/null @@ -1,68 +0,0 @@ -% Copyright (C) 2001 Michel Juillard -% -function dr=resol1(ys,algo,linear,iorder) - -global M_ options_ oo_ bayestopt_ -global it_ means_ stderrs_ dr1_test_ - -dr1_test_ = zeros(2,1); - -if linear == 1 - iorder =1; -end - -xlen = M_.maximum_lead + M_.maximum_lag + 1; -klen = M_.maximum_lag + M_.maximum_lead + 1; -iyv = M_.lead_lag_incidence'; -iyv = iyv(:); -iyr0 = find(iyv) ; -it_ = M_.maximum_lag + 1 ; - -if M_.exo_nbr == 0 - oo_.exo_steady_state = [] ; -end - -if ~ M_.lead_lag_incidence(M_.maximum_lag+1,:) > 0 - error ('RESOL: Error in model specification: some variables don"t appear as current') ; -end - -%if xlen > 1 -% error (['RESOL: stochastic exogenous variables must appear only at the' ... -% ' current period. Use additional endogenous variables']) ; -%end - -% check if ys is steady state -tempex = oo_.exo_simul; -oo_.exo_simul = repmat(transpose(oo_.exo_steady_state), M_.maximum_lag + ... - M_.maximum_lead+1,1); -fh = str2func([M_.fname '_static']); -if max(abs(feval(fh,ys,oo_.exo_simul))) > options_.dynatol - % dirty trick to call either user function or dynare_solve - [dr.ys, cheik] = feval(bayestopt_.static_solve,[M_.fname '_static'],ys,oo_.exo_simul); - if cheik - dr1_test_(1) = 1; % dynare_solve did not converge to the steady state. - resid = feval([M_.fname '_static'],dr.ys,oo_.exo_simul); - dr1_test_(2) = resid'*resid; % penalty... - disp('dynare_solve is unable to find the steady state.') - return - end -else - dr.ys = ys; -end - -dr.fbias = zeros(M_.endo_nbr,1); -dr = dr11(iorder,dr,0); - -if algo == 1 & iorder > 1 - dr.ys = dynare_solve('dr2',ys,dr); - dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_simul); - dr = dr11(iorder,dr,0); -end -oo_.exo_simul = tempex; -tempex = []; - -% 01/01/2003 MJ added dr_algo == 1 -% 08/24/2001 MJ uses Schmitt-Grohe and Uribe (2001) constant correction -% in dr.ghs2 -% 05/26/2003 MJ added temporary values for oo_.exo_simul -% 01/22/2005 SA check --> cheik (to avoid confusion with the function check.m)