function [competitive_rule, dist,exitflag] = competitive_optimal %[competitive_rule, dist,exitflag] = competitive_optimal Computes the optimal operational interest-rate rule coefficients shown in the first row of table 4 of the paper ``Optimal Inflation Stabilization in a Medium-Scale Macroeconomic Model,'' by Stephanie Schmitt-Grohe and Martin Uribe (November 2005). competitive_rule is a 4x1 vector containing the policy rule coefficients (ALFAPAI, ALFAW, ALFAY, and ALFAR); dist is the negative of the second derivative of welfare with respect to the parameter sigma scaling the standard deviation of shocks (i.e., dist=-gss(1)); exitflag equals 1 if the optimization procedure converged (write help fminsearch for more information on exitflag) %Get symbolic expressions for first- and second-order derivatives of equilibrium conditions [f, fx, fy, fxp, fyp, statevar_cu, controlvar_cu, fypyp, fypy, fypxp, fypx, fyyp, fyy, fyxp, fyx, fxpyp, fxpy, fxpxp, fxpx, fxyp, fxy, fxxp, fxx] = competitive; %Ramsey steady state ramsey_ss; load ramsey_ss RSTAR = r_cu; PAISTAR = pai_cu; %find the steady state of welfare h_cu = stil_cup * hd_cu; muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu; uf_cu = log(c_cu - B * c_ba1/muzstar_cu)*(1-PHI4)+log(1-h_cu)*PHI4; vt_cu = uf_cu/(1- BETTA) ; vt_cup = vt_cu; clear h_cu uf_cu muzstar_cu %make all variables into logs u_cu =log(u_cu); g_cu = log(g_cu); muz_cu = log(muz_cu); muupsilon_cu=log(muupsilon_cu); qq_cu = log(qq_cu); r_cu = log(r_cu); s_cu = log(s_cu); stil_cu = log(stil_cu); hd_cu = log(hd_cu); k_cu = log(k_cu); iv_cu = log(iv_cu); w_cu = log(w_cu); output_cu = log(output_cu); c_cu = log(c_cu); la_cu = log(la_cu); pai_cu = log(pai_cu); f2_cu = log(f2_cu); x2_cu = log(x2_cu); u_cup =log(u_cup); muz_cup = log(muz_cup); muupsilon_cup=log(muupsilon_cup); qq_cup = log(qq_cup); r_cup = log(r_cup); s_cup = log(s_cup); hd_cup = log(hd_cup); stil_cup = log(stil_cup); k_cup = log(k_cup); iv_cup = log(iv_cup); w_cup = log(w_cup); output_cup = log(output_cup); g_cup = log(g_cup); c_cup = log(c_cup); la_cup = log(la_cup); pai_cup = log(pai_cup); f2_cup = log(f2_cup); x2_cup = log(x2_cup); u_ba1 =log(u_ba1); muz_ba1 = log(muz_ba1); muupsilon_ba1 = log(muupsilon_ba1); qq_ba1 = log(qq_ba1); r_ba1 = log(r_ba1); s_ba1 = log(s_ba1); hd_ba1 = log(hd_ba1); stil_ba1 = log(stil_ba1); k_ba1 = log(k_ba1); iv_ba1 = log(iv_ba1); w_ba1 = log(w_ba1); output_ba1 = log(output_ba1); g_ba1 = log(g_ba1); c_ba1 = log(c_ba1); la_ba1 = log(la_ba1); pai_ba1 = log(pai_ba1); f2_ba1 = log(f2_ba1); x2_ba1 = log(x2_ba1); u_ba1p =log(u_ba1p); muz_ba1p = log(muz_ba1p); muupsilon_ba1p = log(muupsilon_ba1p); qq_ba1p = log(qq_ba1p); r_ba1p = log(r_ba1p); s_ba1p = log(s_ba1p); hd_ba1p = log(hd_ba1p); stil_ba1p = log(stil_ba1p); k_ba1p = log(k_ba1p); iv_ba1p = log(iv_ba1p); w_ba1p = log(w_ba1p); output_ba1p = log(output_ba1p); g_ba1p = log(g_ba1p); c_ba1p = log(c_ba1p); la_ba1p = log(la_ba1p); pai_ba1p = log(pai_ba1p); f2_ba1p = log(f2_ba1p); x2_ba1p = log(x2_ba1p); u_fu1 =log(u_fu1); muz_fu1 = log(muz_fu1); muupsilon_fu1 = log(muupsilon_fu1); qq_fu1 = log(qq_fu1); r_fu1 = log(r_fu1); s_fu1 = log(s_fu1); hd_fu1 = log(hd_fu1); stil_fu1 = log(stil_fu1); k_fu1 = log(k_fu1); iv_fu1 = log(iv_fu1); w_fu1 = log(w_fu1); output_fu1 = log(output_fu1); g_fu1 = log(g_fu1); c_fu1 = log(c_fu1); la_fu1 = log(la_fu1); pai_fu1 = log(pai_fu1); f2_fu1 = log(f2_fu1); x2_fu1 = log(x2_fu1); u_fu1p =log(u_fu1p); muz_fu1p = log(muz_fu1p); muupsilon_fu1p = log(muupsilon_fu1p); qq_fu1p = log(qq_fu1p); r_fu1p = log(r_fu1p); s_fu1p = log(s_fu1p); hd_fu1p = log(hd_fu1p); stil_fu1p = log(stil_fu1p); k_fu1p = log(k_fu1p); iv_fu1p = log(iv_fu1p); w_fu1p = log(w_fu1p); output_fu1p = log(output_fu1p); g_fu1p = log(g_fu1p); c_fu1p = log(c_fu1p); la_fu1p = log(la_fu1p); pai_fu1p = log(pai_fu1p); f2_fu1p = log(f2_fu1p); x2_fu1p = log(x2_fu1p); %Allow for second-order derivatives of Ramsey equilibrium conditions approx = 2; %Initial guess for coefficeints of the optimal operational interest-rate rule MONET_POL0 = [1.5 0 0 0 ]; %Give the policy coefficients ALFA_PAI = MONET_POL0(1); ALFA_W = MONET_POL0(2); ALFA_Y = MONET_POL0(3); ALFA_R = MONET_POL0(4); num_eval eta = zeros(length(statevar_cu),3); eta(end-2:end,:) = [STD_EPSG 0 0; 0 STD_EPSMUZ 0; 0 0 STD_EPSMUUPSILON]; save alles.mat eta THETA nfx nfxp nfy nfyp nfypyp nfypy nfypxp nfypx nfyyp nfyy nfyxp nfyx nfxpyp nfxpy nfxpxp nfxpx nfxyp nfxy nfxxp nfxx OPTIONS=optimset('MaxFunEvals',10000,'MaxIter',10000,'TolFun',1e-11,'TolX',1e-11); [competitive_rule, dist,exitflag] = fminsearch(@(x) get_gap(x), MONET_POL0, OPTIONS); %WELFARE.M function welfare_conditional = get_gap(x); ALFA_PAI = x(1); ALFA_W = x(2); ALFA_Y = x(3); ALFA_R = x(4); load alles.mat %Reevaluate the relevant equations. nfy(16,:) = [0 0 0 ALFA_Y 0 0 ALFA_PAI+ALFA_W 0 0 0 -1 0 ALFA_W 0]; nfx(16,:) = [ALFA_R 0 0 0 -ALFA_W 0 0 0 -ALFA_Y 0 ALFA_W+ALFA_Y ALFA_W*THETA/(1-THETA)+ALFA_Y*THETA/(1-THETA)]; [gx,hx] = gx_hx(nfy,nfx,nfyp,nfxp); if gx ~= exp(1) [gxx,hxx] = gxx_hxx_noloop(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx); [gss,hss] = gss_hss(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx,gxx,eta); welfare_conditional = -gss(1) else welfare_conditional=10; end %gx ~= exp(1) %competitive.M %Equilibrium conditions of the model in ``Optimal Inflation Stabilization in a Medium-Scale Macroeconomic Model'' by S. Schmitt-Grohe and Martin Uribe (November 2005) %(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2005 %Calls allsyms.m, gx_hx.m, ramsey_ss.m, anal_deriv.m, num_eval.m, gss_hss.m, %Subfunctions included in this program: competitive.m, gxx_hxx_noloop.m, hessian.m, solab.m function [f, fx, fy, fxp, fyp, statevar_cu, controlvar_cu, fypyp, fypy, fypxp, fypx, fyyp, fyy, fyxp, fyx, fxpyp, fxpy, fxpxp, fxpx, fxyp, fxy, fxxp, fxx] = competitive; %Create symbolic variables allsyms %Competitive Equilibrium conditions for a given interest-rate rule %Define muzstar_cu, muI_cu muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu; muzstar_cup = muupsilon_cup^(THETA/(1-THETA))*muz_cup; muI_cu = muupsilon_cu*muzstar_cu; muI_cup = muupsilon_cup*muzstar_cup; mula_cu = 1/muzstar_cu;%holds only if phi3=1 mula_cup = 1/muzstar_cup; %Production Function pf_cu = (u_cu*k_cu/muI_cu)^THETA * hd_cu^(1-THETA) - PSSI; %Capital utilization cost a_cu = GAMA1 * (u_cu-1) + GAMA2/2 * (u_cu-1)^2; a_cup = GAMA1 * (u_cup-1) + GAMA2/2 * (u_cup-1)^2; %Capital adjustment costs qiv_cu = iv_cu * (1 - KAPA /2 * (iv_cu/iv_ba1* muI_cu - MUI)^2); d1qiv_cu = qiv_cu/iv_cu - KAPA*( muI_cu * iv_cu/iv_ba1-MUI )*(muI_cu * iv_cu / iv_ba1); d2qiv_cup = (muI_cup * iv_cup/iv_cu)^2 * KAPA * (muI_cup*iv_cup / iv_cu -MUI); rk_cu = diff(a_cu,'u_cu'); rk_cup = diff(a_cup,'u_cup'); v_cu = sqrt(PHI2/PHI1 + 1/PHI1 * (r_cu-1) / r_cu); ell_cu = PHI1 * v_cu + PHI2 / v_cu - 2 * sqrt(PHI1*PHI2); %transaction cost mh_cu = c_cu / v_cu; mc_cu = rk_cu / THETA / (u_cu*k_cu /hd_cu / muI_cu )^(THETA-1); ptil_cu = ((1 - ALFA * pai_cu ^(ETA-1) * pai_ba1 ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA)); ptil_cup = ((1 - ALFA * pai_cup^(ETA-1) * pai_ba1p^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA)); x1_cu = (ETA -1 ) / ETA * x2_cu; x1_cup = (ETA -1 ) / ETA * x2_cup; wtil_cu = ((w_cu ^(1-ETATIL) - ALFATIL * (w_ba1/muzstar_cu)^(1-ETATIL) * ((MUZSTAR*pai_ba1)^CHITIL / pai_cu) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL)); wtil_cup = ((w_cup^(1-ETATIL) - ALFATIL * (w_ba1p/muzstar_cup) ^(1-ETATIL) * ((MUZSTAR*pai_ba1p)^CHITIL / pai_cup) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL)); h_cu = stil_cup * hd_cu; h_cup = stil_fu1p * hd_cup; f1_cu = f2_cu; f1_cup = f2_cup; %Competitive equilibrium Conditions e1_cu = -k_cup + (1-DELTA) * k_cu/muI_cu + qiv_cu; e2_cu = -la_cu * (1+2*PHI1*v_cu-2*(PHI1*PHI2)^(1/2)) + (1-PHI4) * (c_cu - B * c_ba1/muzstar_cu)^((1-PHI4)*(1-PHI3)-1) * (1-h_cu)^(PHI4*(1-PHI3)) - B * BETTA * (1-PHI4) * (muzstar_cup*c_cup - B * c_cu)^((1-PHI4)*(1-PHI3)-1) * (1-h_cup)^(PHI4*(1-PHI3)); e3_cu = -la_cu * qq_cu + BETTA * mula_cup/muupsilon_cup * la_cup * ((u_cup *rk_cup - a_cup) + qq_cup * (1-DELTA )); e4_cu = -la_cu + la_cu * qq_cu * d1qiv_cu + BETTA * mula_cup / muupsilon_cup * la_cup * qq_cup * d2qiv_cup; e5_cu = -f1_cu + (ETATIL-1) / ETATIL * wtil_cu * la_cu * (w_cu/wtil_cu)^ETATIL * hd_cu + ALFATIL * BETTA * (pai_cup / (MUZSTAR*pai_cu)^CHITIL)^(ETATIL-1) * (muzstar_cup*wtil_cup/wtil_cu)^(ETATIL-1) * mula_cup * muzstar_cup * f1_cup; e6_cu = -f2_cu + PHI4 * (c_cu-B*c_ba1/muzstar_cu)^((1-PHI3)*(1-PHI4)) * (1-h_cu)^(PHI4*(1-PHI3)-1) * (w_cu/wtil_cu)^ETATIL * hd_cu + ALFATIL * BETTA * (pai_cup / (MUZSTAR * pai_cu)^CHITIL)^ETATIL * mula_cup * muzstar_cup * f2_cup* (muzstar_cup* wtil_cup/wtil_cu)^ETATIL; e7_cu = -la_cu + BETTA * r_cu * mula_cup * la_cup / pai_cup; e8_cu = -x1_cu + output_cu * mc_cu * ptil_cu^(-ETA-1) + ALFA * BETTA * mula_cup * la_cup / la_cu * (ptil_cu / ptil_cup)^(-ETA-1) * (pai_cu^CHI / pai_cup)^(-ETA) * muzstar_cup * x1_cup; e9_cu = -x2_cu + output_cu * ptil_cu^(-ETA) + ALFA * BETTA * mula_cup * la_cup / la_cu * (pai_cu^CHI / pai_cup )^(1-ETA) * (ptil_cu / ptil_cup)^(-ETA) * muzstar_cup * x2_cup; e10_cu = -pf_cu + output_cu * s_cup; e11_cu = - s_cup + (1 - ALFA) * ptil_cu^(-ETA) + ALFA * (pai_cu/ pai_ba1^CHI)^ETA * s_cu; e12_cu = - stil_cup + (1-ALFATIL) * (wtil_cu/w_cu)^(-ETATIL) + ALFATIL * (w_ba1 / w_cu/ muzstar_cu ) ^(-ETATIL) * (pai_cu / (MUZSTAR * pai_ba1)^CHITIL)^ETATIL * stil_cu; %Investment e14_cu = -iv_cu + output_cu - c_cu * (1+ell_cu) - g_cu - a_cu * k_cu / muI_cu; %real wage e15_cu = -w_cu + mc_cu * (1-THETA) * (u_cu * k_cu / muI_cu / hd_cu )^(THETA) / (1 + NU * (r_cu -1)/r_cu ); %Utility function uf_cu = log(c_cu - B * c_ba1 / muzstar_cu)*(1-PHI4)+log(1-h_cu)*PHI4; %Value function e16_cu = -vt_cu + uf_cu + BETTA * vt_cup; %Monetary Policy Rule e17_cu = -log(r_cu/RSTAR) + ALFA_R * log(r_ba1/RSTAR)+ ALFA_PAI * log(pai_cu/PAISTAR) + ALFA_W *( log(w_cu/w_ba1) + log(pai_cu/PAISTAR) + log(muzstar_cu/MUZSTAR)) + ALFA_Y * log(output_cu/output_ba1*muzstar_cu/MUZSTAR); %Transform equilibrium conditions (which is a system of higher-than-first order nonlinear difference equations) into a system of first-order nonlinear difference equations e18_cu = [- r_ba1p + r_cu; - c_ba1p + c_cu; - iv_ba1p + iv_cu; - pai_ba1p + pai_cu; - stil_cup + stil_fu1; - w_ba1p + w_cu; -output_ba1p+output_cu]; %Place all equilibrium conditions into a vector constraints_cu = [e1_cu; e2_cu; e3_cu; e4_cu; e5_cu; -e6_cu; e7_cu; -e8_cu; e9_cu; -e10_cu; -e11_cu; -e12_cu; e14_cu; e15_cu; e16_cu; e17_cu; e18_cu]; %Government Purchases shock exprocg = log(g_cup/G) - RHOG * log(g_cu/G); %Evolution of Neutral Technolgy shock exprocmuz = log(muz_cup/MUZ) - RHOMUZ * log(muz_cu/MUZ); %Evolution of Investment-Specific Technolgy shock exprocmuupsilon = log(muupsilon_cup/MUUPSILON) - RHOMUUPSILON * log(muupsilon_cu/MUUPSILON); %stack exogenous processes into a vector exproc = [exprocg;exprocmuz;exprocmuupsilon]; %Complete set of equilibrium conditions a a system of first-order nonlinear difference equations f = [constraints_cu; exproc]; %state varibles statevar_cu = [r_ba1 c_ba1 iv_ba1 pai_ba1 w_ba1 s_cu stil_cu k_cu output_ba1 g_cu muz_cu muupsilon_cu]; %state variables in t+1 statevar_cup = [r_ba1p c_ba1p iv_ba1p pai_ba1p w_ba1p s_cup stil_cup k_cup output_ba1p g_cup muz_cup muupsilon_cup]; %control variables controlvar_cu = [vt_cu qq_cu hd_cu output_cu c_cu la_cu pai_cu u_cu f2_cu x2_cu r_cu iv_cu w_cu stil_fu1 ]; %control variables in t+1 controlvar_cup = [vt_cup qq_cup hd_cup output_cup c_cup la_cup pai_cup u_cup f2_cup x2_cup r_cup iv_cup w_cup stil_fu1p]; %Make f a function of the logarithm of the state and control vector %controls in logs lc = find(controlvar_cu~='vt_cu'); %variables to substitute from levels to logs vs=[statevar_cu, controlvar_cu(lc), statevar_cup, controlvar_cup(lc)]; f = subs(f, vs, exp(vs)); %Allow for a second-order approximation of the equilibrium approx = 2; %Compute analytical derivatives of f [fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,statevar_cu,controlvar_cu,statevar_cup,controlvar_cup,approx); function [gxx,hxx] = gxx_hxx_noloop(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx) %This program finds the second drivatives of the funtions g(x) and h(x), denoted by gxx and hxx respectively. %It solves the following system: % A + B g_{xx} h_x + C g_{xx} + D h_{xx} =0 nfypyp=hessian(nfypyp); nfypy=hessian(nfypy); nfypxp=hessian(nfypxp); nfypx=hessian(nfypx); nfyyp=hessian(nfyyp); nfyy =hessian(nfyy); nfyxp=hessian(nfyxp); nfyx =hessian(nfyx); nfxpyp=hessian(nfxpyp); nfxpy=hessian(nfxpy); nfxpxp=hessian(nfxpxp); nfxpx=hessian(nfxpx); nfxyp=hessian(nfxyp); nfxy=hessian(nfxy); nfxxp=hessian(nfxxp); nfxx=hessian(nfxx); nx=size(nfx,2); ny = size(nfy,2); n = nx + ny; A = kron(eye(n),hx'*gx')*(nfypyp*gx*hx + nfypxp*hx + nfypy*gx + nfypx) + ... kron(eye(n),gx') *(nfyyp *gx*hx + nfyxp *hx + nfyy *gx + nfyx ) + ... kron(eye(n),hx') *(nfxpyp*gx*hx + nfxpxp*hx + nfxpy*gx + nfxpx) + ... (nfxyp *gx*hx + nfxxp *hx + nfxy *gx + nfxx); B = kron(nfyp, hx') ; C = kron(nfy, eye(nx)); D = kron(nfyp*gx, eye(nx))+ kron(nfxp, eye(nx)); Qq = -[ kron(hx',B)+kron(eye(nx),C), kron(eye(nx), D)] \ (A(:)); gxx = permute(reshape(Qq(1:nx^2*ny),nx,ny,nx),[2 1 3]); hxx = permute(reshape(Qq(nx^2*ny+1:end),nx,nx,nx),[2 1 3]); function x = hessian(x) [n1,n2,n3]=size(x); x=reshape(permute(x,[ 2 1 3]), n1*n2, n3); % % Function: solab % % Purpose: Solves for the recursive representation of the stable solution to a system % of linear difference equations. % % Inputs: Two square matrices a and b and a natural number nk % % a and b are the coefficient matrices of the difference equation % % a*x(t+1) = b*x(t) % % where x(t) is arranged so that the state variables come first, and % % nk is the number of state variables. % % Outputs: the decision rule f and the law of motion p. If we write % % x(t) = [k(t);u(t)] where k(t) contains precisely the state variables, then % % u(t) = f*k(t) and % % k(t+1) = p*k(t). % % Calls: qzdiv % %This program is a slight modificaiton of Paul Klein's origianl solab.m function [f,p] = solab(a,b,nk); [s,t,q,z] = qz(a,b); % upper triangular factorization of the matrix pencil b-za [s,t,q,z] = qzdiv(1,s,t,q,z); % reordering of generalized eigenvalues in ascending order z21 = z(nk+1:end,1:nk); z11 = z(1:nk,1:nk); if rank(z11)abs(s(nk,nk)) | abs(t(nk+1,nk+1))abs(s(nk,nk)) | abs(t(nk+1,nk+1))