From eb1a26962e9f9d0dd23268ddf2b0f191c094acfa Mon Sep 17 00:00:00 2001 From: assia Date: Mon, 3 Mar 2008 09:30:06 +0000 Subject: [PATCH] isolated functions deleted git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1727 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/dr11.m | 255 ---------------------------- matlab/fbeta.m | 6 - matlab/fgamma.m | 4 - matlab/figamm.m | 2 - matlab/lpdfbeta.m | 24 --- matlab/posterior_density_estimate.m | 178 ------------------- 6 files changed, 469 deletions(-) delete mode 100644 matlab/dr11.m delete mode 100644 matlab/fbeta.m delete mode 100644 matlab/fgamma.m delete mode 100644 matlab/figamm.m delete mode 100644 matlab/lpdfbeta.m delete mode 100644 matlab/posterior_density_estimate.m diff --git a/matlab/dr11.m b/matlab/dr11.m deleted file mode 100644 index 18bad98e9..000000000 --- a/matlab/dr11.m +++ /dev/null @@ -1,255 +0,0 @@ -% Copyright (C) 2001 Michel Juillard -% -function dr=dr11(iorder,dr,cheik) - -global M_ options_ oo_ -global it_ stdexo_ means_ dr1_test_ bayestopt_ - -% hack for Bayes -global dr1_test_ bayestopt_ - -options_ = set_default_option(options_,'loglinear',0); - -xlen = M_.maximum_lead + M_.maximum_lag + 1; -klen = M_.maximum_lag + M_.maximum_lead + 1; -iyv = transpose(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 ('Error in model specification: some variables don"t appear as current') ; -end - -if ~cheik -% if xlen > 1 -% error (['SS: stochastic exogenous variables must appear only at the' ... -% ' current period. Use additional endogenous variables']) ; -% end -end - -if M_.maximum_lead > 1 & iorder > 1 - error (['Models with leads on more than one period can only be solved' ... - ' at order 1']) -end - -dr=set_state_space(dr); -kstate = dr.kstate; -kad = dr.kad; -kae = dr.kae; -nstatic = dr.nstatic; -nfwrd = dr.nfwrd; -npred = dr.npred; -nboth = dr.nboth; -order_var = dr.order_var; -nd = size(kstate,1); - -sdyn = M_.endo_nbr - nstatic; - - -tempex = oo_.exo_simul; - -it_ = M_.maximum_lag + 1; -z = repmat(dr.ys,1,klen); -z = z(iyr0) ; -%M_.jacobia=real(diffext('ff1_',[z; oo_.exo_steady_state])) ; -%M_.jacobia=real(jacob_a('ff1_',[z; oo_.exo_steady_state])) ; -[junk,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul); -oo_.exo_simul = tempex ; -tempex = []; - -nz = size(z,1); -k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_lag+1),:); -b = M_.jacobia(:,M_.lead_lag_incidence(M_.maximum_lag+1,order_var)); -a = b\M_.jacobia(:,nonzeros(k1')); -if any(isinf(a(:))) - dr1_test_(1) = 5; - dr1_test_(2) = bayestopt_.penalty; -end -if M_.exo_nbr - fu = b\M_.jacobia(:,nz+1:end); -end - -if M_.maximum_lead == 0 & M_.maximum_lag == 1; % backward model with one lag - dr.ghx = -a; - dr.ghu = -fu; - return; -elseif M_.maximum_lead == 0 & M_.maximum_lag > 1 % backward model with lags on more than - % one period - e = zeros(endo_nbr,nd); - k = find(kstate(:,2) <= M_.maximum_lag+1 & kstate(:,4)); - e(:,k) = -a(:,kstate(k,4)) ; - dr.ghx = e; - dr.ghu = -fu; -end - -% buildind D and E -d = zeros(nd,nd) ; -e = d ; - -k = find(kstate(:,2) >= M_.maximum_lag+2 & kstate(:,3)); -d(1:sdyn,k) = a(nstatic+1:end,kstate(k,3)) ; -k1 = find(kstate(:,2) == M_.maximum_lag+2); -a1 = eye(sdyn); -e(1:sdyn,k1) = -a1(:,kstate(k1,1)-nstatic); -k = find(kstate(:,2) <= M_.maximum_lag+1 & kstate(:,4)); -e(1:sdyn,k) = -a(nstatic+1:end,kstate(k,4)) ; -k2 = find(kstate(:,2) == M_.maximum_lag+1); -k2 = k2(~ismember(kstate(k2,1),kstate(k1,1))); -d(1:sdyn,k2) = a1(:,kstate(k2,1)-nstatic); - -if ~isempty(kad) - for j = 1:size(kad,1) - d(sdyn+j,kad(j)) = 1 ; - e(sdyn+j,kae(j)) = 1 ; - end -end -options_ = set_default_option(options_,'qz_criterium',1.000001); - -if ~exist('mjdgges') - % using Chris Sim's routines - use_qzdiv = 1; - [ss,tt,qq,w] = qz(e,d); - [tt,ss,qq,w] = qzdiv(options_.qz_criterium,tt,ss,qq,w); - ss1=diag(ss); - tt1=diag(tt); - warning_state = warning; - warning off; - oo_.eigenvalues = ss1./tt1 ; - warning warning_state; - nba = nnz(abs(eigval) > options_.qz_criterium); -else - use_qzdiv = 0; - [ss,tt,w,sdim,oo_.eigenvalues,info] = mjdgges(e,d,options_.qz_criterium); - if info & info ~= nd+2; - error(['ERROR' info ' in MJDGGES.DLL']); - end - nba = nd-sdim; -end - -nyf = sum(kstate(:,2) > M_.maximum_lag+1); - -if cheik - dr.rank = rank(w(1:nyf,nd-nyf+1:end)); - % dr.eigval = oo_.eigenvalues; - return -end - -eigenvalues = sort(oo_.eigenvalues); - -if nba > nyf; -% disp('Instability !'); - dr1_test_(1) = 3; %% More eigenvalues superior to unity than forward variables ==> instability. - dr1_test_(2) = (abs(eigenvalues(nd-nba+1:nd-nyf))-1-1e-5)'*... - (abs(eigenvalues(nd-nba+1:nd-nyf))-1-1e-5);% Distance to Blanchard-Khan conditions (penalty) - return -elseif nba < nyf; -% disp('Indeterminacy !'); - dr1_test_(1) = 2; %% ==> Indeterminacy. - dr1_test_(2) = (abs(eigenvalues(nd-nyf+1:nd-nba))-1-1e-5)'*... - (abs(eigenvalues(nd-nyf+1:nd-nba))-1-1e-5);% Distance to Blanchard-Khan conditions (penality) - %% warning('DR1: Blanchard-Kahn conditions are not satisfied. Run CHEIK to learn more!'); - return -end - -np = nd - nyf; -n2 = np + 1; -n3 = nyf; -n4 = n3 + 1; -% derivatives with respect to dynamic state variables -% forward variables - -if condest(w(1:n3,n2:nd)) > 1e9 -% disp('Indeterminacy !!'); - dr1_test_(1) = 2; - dr1_test_(2) = 1; - return -end - -warning_state = warning; -lastwarn(''); -warning off; -gx = -w(1:n3,n2:nd)'\w(n4:nd,n2:nd)'; - -if length(lastwarn) > 0; -% disp('Indeterminacy !!'); - dr1_test_(1) = 2; - dr1_test_(2) = 1; - warning(warning_state); - return -end - -% 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); - -lastwarn(''); -if length(lastwarn) > 0; -% disp('Singularity problem in dr11.m'); - dr1_test_(1) = 2; - dr1_test_(2) = 1; - warning(warning_state); - return -end - -k1 = find(kstate(n4:nd,2) == M_.maximum_lag+1); -k2 = find(kstate(1:n3,2) == M_.maximum_lag+2); -dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)]; - -%lead variables actually present in the model -j3 = nonzeros(kstate(:,3)); -j4 = find(kstate(:,3)); -% derivatives with respect to exogenous variables -if M_.exo_nbr - a1 = eye(M_.endo_nbr); - aa1 = []; - if nstatic > 0 - aa1 = a1(:,1:nstatic); - end - dr.ghu = -[aa1 a(:,j3)*gx(j4,1:npred)+a1(:,nstatic+1:nstatic+ ... - npred) a1(:,nstatic+npred+1:end)]\fu; - - - lastwarn(''); - if length(lastwarn) > 0; -% disp('Singularity problem in dr11.m'); - dr1_test_(1) = 2; - dr1_test_(2) = 1; - return - end -end -warning(warning_state); - -% static variables -if nstatic > 0 - 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 - -if options_.loglinear == 1 - k = find(dr.kstate(:,2) <= M_.maximum_lag+1); - klag = dr.kstate(k,[1 2]); - k1 = dr.order_var; - - dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ... - repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1); - dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu; -end - -% necessary when using Sims' routines -if use_qzdiv - gx = real(gx); - hx = real(hx); - dr.ghx = real(dr.ghx); - dr.ghu = real(dr.ghu); -end - - diff --git a/matlab/fbeta.m b/matlab/fbeta.m deleted file mode 100644 index 75c8d055f..000000000 --- a/matlab/fbeta.m +++ /dev/null @@ -1,6 +0,0 @@ -function e = fbeta(p2,p,p1,perc) -% must restrict p2 such that a>0 and b>0 .... - a = (1-p1)*p1^2/p2^2 - p1; - b = a*(1/p1 - 1); - - e = p - pbeta(perc,a,b); \ No newline at end of file diff --git a/matlab/fgamma.m b/matlab/fgamma.m deleted file mode 100644 index 209771148..000000000 --- a/matlab/fgamma.m +++ /dev/null @@ -1,4 +0,0 @@ -function e = fgamma(p2,p,p1,perc) - b = p2^2/p1; - a = p1/b; - e = p - pgamma(perc,a,b); \ No newline at end of file diff --git a/matlab/figamm.m b/matlab/figamm.m deleted file mode 100644 index 37d339748..000000000 --- a/matlab/figamm.m +++ /dev/null @@ -1,2 +0,0 @@ -function e = figamm(p2,p,p1,perc) - e = p - pgamma(1/perc,p2/2,2/p1); \ No newline at end of file diff --git a/matlab/lpdfbeta.m b/matlab/lpdfbeta.m deleted file mode 100644 index 13917bbd6..000000000 --- a/matlab/lpdfbeta.m +++ /dev/null @@ -1,24 +0,0 @@ -function ldens = lpdfbeta(x,a,b); - -% function ldens = lpdfbeta(x,a,b) -% log Beta PDF -% -% INPUTS -% x: density evatuated at x -% a: Beta distribution paramerer -% b: Beta distribution paramerer -% -% OUTPUTS -% ldens: the log Beta PDF -% -% SPECIAL REQUIREMENTS -% none -% -% part of DYNARE, copyright Dynare Team (2003-2008) -% Gnu Public License. - - ldens = gammaln(a+b) - gammaln(a) - gammaln(b) + (a-1)*log(x) + (b-1)*log(1-x); - -% 10/11/03 MJ adapted from a GAUSS version by F. Schorfheide, translated -% to Matlab by R. Wouters. -% use Matlab gammaln instead of lngam \ No newline at end of file diff --git a/matlab/posterior_density_estimate.m b/matlab/posterior_density_estimate.m deleted file mode 100644 index 8ac63317e..000000000 --- a/matlab/posterior_density_estimate.m +++ /dev/null @@ -1,178 +0,0 @@ -function [abscissa,f,h] = posterior_density_estimate(data,number_of_grid_points,bandwidth,kernel_function) -%% -%% This function aims at estimating posterior univariate densities from realisations of a Metropolis-Hastings -%% algorithm. A kernel density estimator is used (see Silverman [1986]) and the main task of this function is -%% to obtain an optimal bandwidth parameter. -%% -%% * Silverman [1986], "Density estimation for statistics and data analysis". -%% * M. Skold and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm". -%% -%% The last section is adapted from Anders Holtsberg's matlab toolbox (stixbox). -%% -%% stephane.adjemian@cepremap.cnrs.fr [01/16/2004]. - -if size(data,2) > 1 & size(data,1) == 1; - data = data'; -elseif size(data,2)>1 & size(data,1)>1; - error('density_estimate: data must be a one dimensional array.'); -end; -test = log(number_of_grid_points)/log(2); -if ( abs(test-round(test)) > 10^(-12)); - error('The number of grid points must be a power of 2.'); -end; - -n = size(data,1); - - -%% KERNEL SPECIFICATION... - -if strcmp(kernel_function,'gaussian'); - k = inline('inv(sqrt(2*pi))*exp(-0.5*x.^2)'); - k2 = inline('inv(sqrt(2*pi))*(-exp(-0.5*x.^2)+(x.^2).*exp(-0.5*x.^2))'); % second derivate of the gaussian kernel - k4 = inline('inv(sqrt(2*pi))*(3*exp(-0.5*x.^2)-6*(x.^2).*exp(-0.5*x.^2)+(x.^4).*exp(-0.5*x.^2))'); % fourth derivate... - k6 = inline('inv(sqrt(2*pi))*(-15*exp(-0.5*x.^2)+45*(x.^2).*exp(-0.5*x.^2)-15*(x.^4).*exp(-0.5*x.^2)+(x.^6).*exp(-0.5*x.^2))'); % sixth derivate... - mu02 = inv(2*sqrt(pi)); - mu21 = 1; -elseif strcmp(kernel_function,'uniform'); - k = inline('0.5*(abs(x) <= 1)'); - mu02 = 0.5; - mu21 = 1/3; -elseif strcmp(kernel_function,'triangle'); - k = inline('(1-abs(x)).*(abs(x) <= 1)'); - mu02 = 2/3; - mu21 = 1/6; -elseif strcmp(kernel_function,'epanechnikov'); - k = inline('0.75*(1-x.^2).*(abs(x) <= 1)'); - mu02 = 3/5; - mu21 = 1/5; -elseif strcmp(kernel_function,'quartic'); - k = inline('0.9375*((1-x.^2).^2).*(abs(x) <= 1)'); - mu02 = 15/21; - mu21 = 1/7; -elseif strcmp(kernel_function,'triweight'); - k = inline('1.09375*((1-x.^2).^3).*(abs(x) <= 1)'); - k2 = inline('(105/4*(1-x.^2).*x.^2-105/16*(1-x.^2).^2).*(abs(x) <= 1)'); - k4 = inline('(-1575/4*x.^2+315/4).*(abs(x) <= 1)'); - k6 = inline('(-1575/2).*(abs(x) <= 1)'); - mu02 = 350/429; - mu21 = 1/9; -elseif strcmp(kernel_function,'cosinus'); - k = inline('(pi/4)*cos((pi/2)*x).*(abs(x) <= 1)'); - k2 = inline('(-1/16*cos(pi*x/2)*pi^3).*(abs(x) <= 1)'); - k4 = inline('(1/64*cos(pi*x/2)*pi^5).*(abs(x) <= 1)'); - k6 = inline('(-1/256*cos(pi*x/2)*pi^7).*(abs(x) <= 1)'); - mu02 = (pi^2)/16; - mu21 = (pi^2-8)/pi^2; -end; - - -%% OPTIMAL BANDWIDTH PARAMETER.... - -if bandwidth == 0; % Rule of thumb bandwidth parameter (Silverman [1986] corrected by - % Skold and Roberts [2003] for Metropolis-Hastings). - sigma = std(data); - h = 2*sigma*(sqrt(pi)*mu02/(12*(mu21^2)*n))^(1/5); % Silverman's optimal bandwidth parameter. - A = 0; - for i=1:n; - j = i; - while j<= n & data(j,1)==data(i,1); - j = j+1; - end; - A = A + 2*(j-i) - 1; - end; - A = A/n; - h = h*A^(1/5); % correction -elseif bandwidth == -1; % Adaptation of the Sheather and Jones [1991] plug-in estimation of the optimal bandwidth - % parameter for metropolis hastings algorithm. - if strcmp(kernel_function,'uniform') | ... - strcmp(kernel_function,'triangle') | ... - strcmp(kernel_function,'epanechnikov') | ... - strcmp(kernel_function,'quartic'); - error('I can''t compute the optimal bandwidth with this kernel... Try the gaussian, triweight or cosinus kernels.'); - end; - sigma = std(data); - A = 0; - for i=1:n; - j = i; - while j<= n & data(j,1)==data(i,1); - j = j+1; - end; - A = A + 2*(j-i) - 1; - end; - A = A/n; - Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi)); - g3 = abs(2*A*k6(0)/(mu21*Itilda4*n))^(1/9); - Ihat3 = 0; - for i=1:n; - Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3)); - end; - Ihat3 = -Ihat3/((n^2)*g3^7); - g2 = abs(2*A*k4(0)/(mu21*Ihat3*n))^(1/7); - Ihat2 = 0; - for i=1:n; - Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2)); - end; - Ihat2 = Ihat2/((n^2)*g2^5); - h = (A*mu02/(n*Ihat2*mu21^2))^(1/5); % equation (22) in Skold and Roberts [2003] --> h_{MH} -elseif bandwidth == -2; % Bump killing... We construct local bandwith parameters in order to remove - % spurious bumps introduced by long rejecting periods. - if strcmp(kernel_function,'uniform') | ... - strcmp(kernel_function,'triangle') | ... - strcmp(kernel_function,'epanechnikov') | ... - strcmp(kernel_function,'quartic'); - error('I can''t compute the optimal bandwidth with this kernel... Try the gaussian, triweight or cosinus kernels.'); - end; - sigma = std(data); - A = 0; - T = zeros(n,1); - for i=1:n; - j = i; - while j<= n & data(j,1)==data(i,1); - j = j+1; - end; - T(i) = (j-i); - A = A + 2*T(i) - 1; - end; - A = A/n; - Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi)); - g3 = abs(2*A*k6(0)/(mu21*Itilda4*n))^(1/9); - Ihat3 = 0; - for i=1:n; - Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3)); - end; - Ihat3 = -Ihat3/((n^2)*g3^7); - g2 = abs(2*A*k4(0)/(mu21*Ihat3*n))^(1/7); - Ihat2 = 0; - for i=1:n; - Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2)); - end; - Ihat2 = Ihat2/((n^2)*g2^5); - h = ((2*T-1)*mu02/(n*Ihat2*mu21^2)).^(1/5); % Note that h is a column vector (local banwidth parameters). -elseif bandwidth > 0; - h = bandwidth; -else; - error('density_estimate: bandwidth must be positive or equal to 0,-1 or -2.'); -end; - -%% COMPUTE DENSITY ESTIMATE, using the optimal bandwidth parameter. -%% -%% This section is adapted from Anders Holtsberg's matlab toolbox -%% (stixbox --> plotdens.m). - - -a = min(data) - (max(data)-min(data))/3; -b = max(data) + (max(data)-min(data))/3; -abscissa = linspace(a,b,number_of_grid_points)'; -d = abscissa(2)-abscissa(1); -xi = zeros(number_of_grid_points,1); -xa = (data-a)/(b-a)*number_of_grid_points; -for i = 1:n; - indx = floor(xa(i)); - temp = xa(i)-indx; - xi(indx+[1 2]) = xi(indx+[1 2]) + [1-temp,temp]'; -end; -xk = [-number_of_grid_points:number_of_grid_points-1]'*d; -kk = k(xk/h); -kk = kk / (sum(kk)*d*n); -f = ifft(fft(fftshift(kk)).*fft([xi ;zeros(size(xi))])); -f = real(f(1:number_of_grid_points)); \ No newline at end of file