% This is the Matlab program for the Econ Enviro Paper % Solution method following UHLIG undetermined coefficients with % sensitivity to the structure as described in Marimon (2001) % COPYRIGHT by O. MIKHAIL (2005) % Feel free to copy, modify and use at your own risk. % However, you are not allowed to sell this code or otherwise impinge % on its free distribution. clear all; clc; diary off; dos('del econ-enviro-phi2.out'); diaryfile = ('econ-enviro-phi2.out'); diary(diaryfile); TT = clock; fprintf(1,'\n'); fprintf(1,'DATE\n'); fprintf(1,'month/day/year %d/%d/%d \n',TT(2),TT(3),TT(1)); fprintf(1,'TIME\n'); fprintf(1,'Hour:minutes:secondes %5.2f:%5.2f:%5.2f \n',TT(4),TT(5),TT(6)); dographs = 1; dolatex = 0; format compact; disp('******************************************************'); disp('******************************************************'); disp('ECON ENVIRO PAPER'); disp(''); disp('Written by O. Mikhail'); disp('January 2005'); disp('Solving using UHLIG with Sensitivity to the Structure'); disp('******************************************************'); disp('******************************************************'); HORIZON = 28; resp_all = zeros(HORIZON,24); ii = 1; jj = 8; phivec = [0.2 0.5 0.8] [notimportant,numberpoints]=size(phivec); for i=1:1:numberpoints phi = phivec(i); disp('******************************************************'); disp('******************************************************'); disp('*********** LOOP OVER PARAMTER PHI *****************'); disp('******************************************************'); disp(sprintf('*********** PHI = %3.3f *************************',phi)) disp('******************************************************'); disp('******************************************************'); disp('******************************************************'); syms cbar ebar mbar lambdabar ybar ibar syms alpha zbar phi beta Abar deltaK deltaE STST = solve('ybar - Abar*kbar^alpha',... 'ibar - deltaK*kbar',... '(1/cbar)-lambdabar-(lambdabar*(phi/zbar))',... '-lambdabar*(1/zbar)+(beta*lambdabar*((1-deltaE)/zbar))+(beta/ebar)',... 'ybar - mbar - cbar - ibar',... 'zbar*mbar - (deltaE*ebar)+(phi*cbar)',... 'kbar = ((1/beta - 1 + deltaK )/(alpha*Abar))^(1/(alpha-1))',... 'cbar','ebar','mbar','lambdabar','ybar','ibar','kbar') cbarf = simple(STST.cbar); ebarf = simple(STST.ebar); mbarf = simple(STST.mbar); lbarf = simple(STST.lambdabar); ybarf = simple(STST.ybar); ibarf = simple(STST.ibar); kbarf = simple(STST.kbar); disp('Hit Enter - Steady State is Solved'); pause; % Input the parameter values beta = 0.99; % beta for quarterly interest rate = 1 percent deltaK = 0.06; % depreciation rate for Capital K deltaEvar = 0.2; % depreciation rate for Enviro Quality E alpha = 0.35; % Capital Share in production zbar = 1.1; % Normalization for the Shock rho = 0.9; % first order correlation for the shock t1 = 1000; % Length to simulate for the time series sigma = 0.001; % The Standard Deviation of the lognormal Shock dum = 1; % dum equals to 0 if you want impulse responses % dum equals to 1 if you want to simulate the series phi = phivec(i); deltaE = deltaEvar; Abar = 1.1; deltaE = deltaEvar; phi = phivec(i); css = subs(STST.cbar); ess = subs(STST.ebar); mss = subs(STST.mbar); lss = subs(STST.lambdabar); yss = subs(STST.ybar); iss = subs(STST.ibar); kss = subs(STST.kbar); disp(' ') disp(' The Parameters Used') fprintf(1,' beta %f \n',beta) fprintf(1,' deltaK %f \n',deltaK) fprintf(1,' deltaE %f \n',deltaE) fprintf(1,' alpha %f \n',alpha) fprintf(1,' rho %f \n',rho) fprintf(1,' phi %f \n',phi) fprintf(1,' rho %f \n',rho) fprintf(1,' Capital Output ratio %f \n',kss/yss) disp(' ') disp(' The Steady State Values ') fprintf(1,' Consumption css %f \n',css) fprintf(1,' Capital kss %f \n',kss) fprintf(1,' Enviro Quality ess %f \n',ess) fprintf(1,' Output yss %f \n',yss) fprintf(1,' Investment iss %f \n',iss) fprintf(1,' Enviro Invest mss %f \n',mss) fprintf(1,' Lambda lss %f \n',lss) disp(' ') disp(' Hit Enter to Continue ') disp(' ') pause; VARNAMES = ['capital ', % Capital 'enviro ', % Environment 'EnvInvest ', % Invest in Environment 'CapInvest ', % Invest in Capital 'cons ', % Consumption 'lambda1 ', % Lagrange Multiplier 'output ', % GDP 'MShock ']; % Shock disp(' ') disp(' DETERMINISTIC EQUATIONS ') % k(t) e(t) m(t) i(t) AA = [ 0 0 0 0 alpha 0 0 0 0 0 mss iss] % k(t-1) e(t-1) m(t-1) i(t-1) BB = [ 0 0 0 0 0 0 0 0 0 0 0 0 ] % c(t) l(t) y(t) CC = [css^(-1) (lss+(lss*zbar^(-1)*phi)) 0 0 0 -1 css 0 -yss ] % z(t) DD = [-(lss*zbar^(-1)*phi) 0 0] disp(' END of DETERMINISTIC EQUATIONS ') disp(' ') pause disp(' ') disp(' EXPECTATIONAL EQUATIONS ') disp(' ') pkk = beta*alpha*Abar*lss*(alpha-1)*kss^(alpha-1); % k(t+1), e(t+1) m(t+1) i(t+1) FF=[pkk 0 0 0 0 -beta*ess^(-1) 0 0 0 ess 0 0 kss 0 0 0] % k(t), e(t) m(t) i(t) GG = [ 0 0 0 0 0 0 0 0 0 -(1-deltaE)*ess -zbar*mss 0 -(1-deltaK)*kss 0 0 -iss] % k(t-1), e(t-1) m(t-1) i(t-1) HH = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] okk = (beta*lss*(1-deltaK))+(beta*alpha*Abar*lss*kss^(alpha-1)); % c(t+1) l(t+1) y(t+1) JJ = [ 0 okk 0 0 beta*(1-deltaE)*lss*zbar^(-1) 0 0 0 0 0 0 0] % c(t) l(t) y(t) KK = [ 0 -lss 0 0 -lss*zbar^(-1) 0 phi*css 0 0 0 0 0] % For z(t+1) LL = [ 0 -beta*(1-deltaE)*lss*zbar^(-1) 0 0] % For z(t) MM = [ 0 lss*zbar^(-1) -zbar*mss 0] pause NN = [rho] sigma_eps = sigma; Sigma = [sigma]; [l_equ,m_states] = size(AA); [l_equ,n_endog ] = size(CC); [l_equ,k_exog ] = size(DD); DISPLAY_IMMEDIATELY = 1; TOL = .000001; PERIOD = 4; GNP_INDEX = 7; DISPLAY_ROOTS = 1; MANUAL_ROOTS = 1; Xi_manual = [6]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% message = ' '; warnings = []; options; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k_exog = min(size(NN)); CC_plus = pinv(CC); CC_0 = (null(CC'))'; Psi_mat = [ zeros(l_equ-n_endog,m_states) FF - JJ*CC_plus*AA ]; Gamma_mat = [ CC_0 * AA JJ*CC_plus*BB - GG + KK*CC_plus*AA ]; Theta_mat = [ CC_0 * BB KK*CC_plus*BB - HH ]; Xi_mat = [ Gamma_mat, Theta_mat eye(m_states), zeros(m_states) ]; Delta_mat = [ Psi_mat, zeros(m_states) zeros(m_states), eye(m_states) ]; [Xi_eigvec,Xi_eigval] = eig(Xi_mat,Delta_mat); [Xi_sortabs,Xi_sortindex] = sort(abs(diag(Xi_eigval))); Xi_sortvec = Xi_eigvec(1:2*m_states,Xi_sortindex); Xi_sortval = diag(Xi_eigval(Xi_sortindex,Xi_sortindex)); Xi_select = Xi_manual; Lambda_mat = diag(Xi_sortval(Xi_select)); Omega_mat = [Xi_sortvec((m_states+1):(2*m_states),Xi_select)]; PP = Omega_mat*Lambda_mat*pinv(Omega_mat); PP_imag = imag(PP); PP = real(PP); RR = - CC_plus*(AA*PP+BB); VV = [ kron(eye(k_exog),AA), kron(eye(k_exog),CC) kron(NN',FF)+kron(eye(k_exog),(FF*PP+JJ*RR+GG)), ... kron(NN',JJ)+kron(eye(k_exog),KK) ]; LLNN_plus_MM = LL*NN + MM; QQSS_vec = - VV \ [ DD(:) LLNN_plus_MM(:) ]; QQ = reshape(QQSS_vec(1:m_states*k_exog),m_states,k_exog); SS = reshape(QQSS_vec((m_states*k_exog+1):((m_states+n_endog)*k_exog)),n_endog,k_exog); WW = [ eye(m_states) , zeros(m_states,k_exog) RR*pinv(PP) , (SS-RR*pinv(PP)*QQ) zeros(k_exog,m_states), eye(k_exog) ]; [m_states,k_exog] = size(QQ); [n_endog,k_exog] = size(SS); disp('Exogenous states z(t):'); disp(VARNAMES((m_states+n_endog+1):(m_states+n_endog+k_exog),:)); disp(' '); disp('Endogenous states x(t):'); disp(VARNAMES(1:m_states,:)); disp(' '); if DISPLAY_ROOTS, disp('All the roots are:'); disp(' root abs(root) '); disp([diag(Xi_eigval(Xi_sortindex,Xi_sortindex)),... abs(diag(Xi_eigval(Xi_sortindex,Xi_sortindex)))] ); disp('The chosen roots are:'); disp(' root abs(root) '); disp([diag(Lambda_mat),abs(diag(Lambda_mat))]); disp(' '); end; disp('PP: Recursive equilibrium law of motion for x(t) on x(t-1):'); disp(PP); disp('QQ: Recursive equilibrium law of motion for x(t) on z(t):'); disp(QQ); disp(' '); disp('Other endogenous variables y(t):'); disp(VARNAMES((m_states+1):(m_states+n_endog),:)); disp(' '); disp('RR: Recursive equilibrium law of motion for y(t) on x(t-1):'); disp(RR); disp('SS: Recursive equilibrium law of motion for y(t) on z(t):'); disp(SS); disp(' '); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% do_it; resp_all(:,ii:jj) = Resp_mat'; nvar = m_states+n_endog+k_exog; ii = ii + (nvar*k_exog); jj = jj + (nvar*k_exog); end; dos('del phi.txt'); save phi.txt resp_all -ASCII -TABS zzz = [ zeros(1,24) ; resp_all ]; resp_all = zzz; impcapital = [resp_all(:,1) resp_all(:,9) resp_all(:,17)]; impenv = [resp_all(:,2) resp_all(:,10) resp_all(:,18)]; impenvm = [resp_all(:,3) resp_all(:,11) resp_all(:,19)]; impinvest = [resp_all(:,4) resp_all(:,12) resp_all(:,20)]; impcons = [resp_all(:,5) resp_all(:,13) resp_all(:,21)]; impoutput = [resp_all(:,7) resp_all(:,15) resp_all(:,23)]; plot(impcons); diary off;