%CEE_RUN.M
%This program computes the second-order approximation to the policy function, that is,   arrays hx, gx, hxx, gxx, hss, gss, and ETAMATRIX, as defined in Schmitt-Grohe and Uribe (JEDC, 2004) for the economy  described in ``Optimal Operational Monetary Policy in the Christiano-Eichenbaum-Evans Model of the U.S. Business Cycle,'' by S.\ Schmitt-Grohe and M.\ Uribe (2004). The interest-rate rule considered is the one Schmitt-Grohe and Uribe (2004) identify as the optimal operational rule. 
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%Date August 18, 2004

function conditional_welfare = cee_run

%Policy parameters
ALFA_PAI = 1.125;
ALFA_Y = -0.0625;
ALFA_R = 0;
%Note: the policy parameter i takes the value -1, indicating that the rule is forward looking. This parameter is given as an imput to CEE_MODEL.M. 

%Read Parameter values
[BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, SIGMAM, KAPA, CHI, CHITIL, GAMA1, GAMA2, PSSI, PHI0, PHI1, PAISTAR, RSTAR,  RHOZ, RHOG, STD_EPSZ, STD_EPSG, Z, QQ, R, PTIL, RK, MUTIL, MC, S, HTIL, HD, STIL, K, H, IV, W, WTIL, OUTPUT, G, MH, C, LA, PAI, U, EXP_VF, za, qq, r, ptil, rk, mutil, mc, s, htil, hd, stil, k, h, iv, w, wtil, output, g, mh, c, la, pai, u, f1, f2, x1, x2, exp_vf, cback, ivback, rback, paiback, outputback, wback, zap, qqp, rp, ptilp, rkp, mutilp, mcp, sp, htilp, hdp, stilp, kp, hp, ivp, wp, wtilp, outputp, gp, mhp, cp, lap, paip, up,  f1p, f2p, x1p, x2p, exp_vfp, cbackp, ivbackp, rbackp, paibackp, outputbackp, wbackp,sig] = cee_ss;

%Value Function at the Non-Stochastic Steady State
VF=exp_vf;

%Derivatives of Equilibrium System 
[fx, fxp, fy, fyp, fypyp, fypy, fypxp, fypx, fyyp, fyy, fyxp, fyx, fxpyp, fxpy, fxpxp, fxpx, fxyp, fxy, fxxp, fxx, f]=cee_model(0,0,1);

%Parameters Defining Innovation Processes 

%Recall that the innovation to the state vector is given by sigma * ETAMATRIX * eps_{t+1}, where sigma is a scalar, ETAMATRIX is a nx*ne matrix and eps_t is a vector of exogenous innovations of size neX1 distributed (0,I);
ne=2; 
ETAMATRIX = zeros(size(fx,2),ne); 
%The gov't purchases and techno shocks take the last 2 elements of the state vector x (see CEE_MODEL.M). Thus,
ETAMATRIX(end-1,1)=STD_EPSG;
ETAMATRIX(end,2)=STD_EPSZ;

approx=2; %Order of approximation of eqm

%Compute numerical derivatives of f
num_eval

%First-Order Approximation
[gx,hx] = gx_hx(nfy,nfx,nfyp,nfxp);

%Second-Order Approximation

[gxx,hxx] = gxx_hxx(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,ETAMATRIX);

conditional_welfare = gss(1)  * sig^2 / 2 + VF

%SUBFUNCTIONS

%CEE_SS.M
%Steady-state of the economy  described in ``Optimal Operational Monetary Policy in the Christiano-Eichenbaum-Evans Model of the U.S. Business Cycle,'' by S.\ Schmitt-Grohe and M.\ Uribe (2004).
function [BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, SIGMAM, KAPA, CHI, CHITIL, GAMA1, GAMA2, PSSI, PHI0, PHI1, PAISTAR, RSTAR,  RHOZ, RHOG, STD_EPSZ, STD_EPSG, Z, QQ, R, PTIL, RK, MUTIL, MC, S, HTIL, HD, STIL, K, H, IV, W, WTIL, OUTPUT, G, MH, C, LA, PAI, U, EXP_VF, za, qq, r, ptil, rk, mutil, mc, s, htil, hd, stil, k, h, iv, w, wtil, output, g, mh, c, la, pai, u, f1, f2, x1, x2, exp_vf, cback, ivback, rback, paiback, outputback, wback, zap, qqp, rp, ptilp, rkp, mutilp, mcp, sp, htilp, hdp, stilp, kp, hp, ivp, wp, wtilp, outputp, gp, mhp, cp, lap, paip, up,  f1p, f2p, x1p, x2p, exp_vfp, cbackp, ivbackp, rbackp, paibackp, outputbackp, wbackp,sig] = cee_ss
%(c) Stephanie Schmitt-Grohe and Martin Uribe, July 1, 2004. 

[BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, SIGMAM, KAPA, CHI, CHITIL, GAMA1, GAMA2, PSSI, PHI0, PHI1, PAISTAR,  RHOZ, RHOG, STD_EPSZ, STD_EPSG, Z, G, sig] = cee_deep;

QQ = 1; %Tobin's q

PAI = PAISTAR; %inflation

RSTAR = PAISTAR / BETTA; %Interest-rate target (gross) 

R = PAI / BETTA; %Nominal interest rate (gross)

PTIL = ((1-ALFA * PAI^((ETA-1)*(1-CHI))) / (1-ALFA)) ^ (1/(1-ETA)); %relative price of optimizing firms

S = (1-ALFA) * PTIL^(-ETA) / (1-ALFA*PAI^((1-CHI) * ETA)); %distortion between output and total production

WTIL_O_W = ((1-ALFATIL * PAI^((CHITIL-1)*(1-ETATIL))) / (1-ALFATIL)) ^ (1/(1-ETATIL)); 

U = sqrt(1 + 2/ GAMA2 * (1/BETTA-1+DELTA-GAMA1)); %Level of capacity utilization

AU = GAMA1*(U-1) + GAMA2/2*(U-1)^2;%OUTPUT LOss from variable capacity utilization 

RK = GAMA1 + GAMA2 * (U-1); %Rental rate of capital

X1_O_X2 = (ETA-1) / ETA;

HTIL_O_HD = WTIL_O_W^(-ETATIL);

STIL = (1-ALFATIL) * WTIL_O_W^(-ETATIL) / (1-ALFATIL * PAI^((1-CHITIL)*ETATIL));%distortion between hours supplied and hours used as factor input

MC = PTIL * (ETA -1 ) * (1- ALFA * BETTA * PAI^(ETA*(1-CHI))) / (ETA * (1-ALFA*BETTA* PAI^((CHI-1)*(1-ETA)))); % marginal cost of producing a unit of final good

H_O_HD = STIL;

MUTIL = (ETATIL / (ETATIL-1)) * ((1-ALFATIL * BETTA * PAI^((1-CHITIL)*(ETATIL-1))) / (1 - ALFATIL * BETTA * PAI^((1-CHITIL)*ETATIL))) * (1/WTIL_O_W); % inverse of wage markup

K_O_HD = (RK / MC / THETA)^(1/(THETA-1)) / U; %aggregate capital factor input

W = MC * (1-THETA) * (U*K_O_HD)^THETA / (1 + NU*(1-1/R));%average real wage rate

WTIL = WTIL_O_W * W; %Wage charged by optimizing unions

%Prepare coefficients to solve quadratic equation in HD

a1 = S * AU * K_O_HD + S * DELTA * K_O_HD - (U*K_O_HD)^THETA;

b1 = S * G + PSSI;

c1 = ((1-B*BETTA) * W * S) / ((1-B) * PHI0 * STIL * MUTIL);

HD = (-b1 - sqrt(b1^2-4*a1*c1)) / (2*a1); %Composite labor input 

HTIL = HTIL_O_HD * HD; %Labor demanded by reoptimizing unions

H = H_O_HD * HD; %Labor supply

K = K_O_HD * HD; %Capital stock

IV = DELTA * K; %investment

C = ((1-B*BETTA) * W ) / ((1-B) * PHI0 * STIL * MUTIL * HD); %Consumption

OUTPUT = C + IV + G + AU * K; %Aggregate absorption

LA = (1-B*BETTA) / C / (1-B); %steady state marginal utility of income

MH = ( LA *(1-BETTA/PAI)/PHI1)^(-1/SIGMAM); %Real balances held by households 

F1 = (LA * WTIL^ETATIL * HTIL / (1 - ALFATIL * BETTA * PAI^((1-CHITIL)*(ETATIL-1)))) / WTIL^ETATIL; 

F2 = (ETATIL-1) / ETATIL * WTIL * F1;

X2 = OUTPUT * PTIL^(-ETA) / (1-ALFA* BETTA * PAI^((CHI-1)*(1-ETA)));

X1 = (ETA-1) / ETA * X2;

UF = log(C*(1-B)) - PHI0/2 * H^2 + PHI1 * MH^(1-SIGMAM)/(1-SIGMAM);

VF = UF/(1-BETTA);

EXP_VF = exp(VF);

za=log(Z); 
qq=log(QQ); 
r=log(R);
ptil=log(PTIL); 
rk=log(RK); 
mutil=log(MUTIL); 
mc=log(MC); 
s=log(S); 
htil=log(HTIL);
hd = log(HD);
stil=log(STIL); 
k=log(K); 
h=log(H); 
iv=log(IV); 
w=log(W); 
wtil=log(WTIL); 
output=log(OUTPUT); 
g=log(G); 
mh=log(MH); 
c=log(C); 
la=log(LA); 
pai=log(PAI); 
u=log(U);
cback=log(C);
ivback=log(IV); 
rback=log(R); 
paiback=log(PAI); 
outputback=log(OUTPUT); 
wback=log(W);
f1=log(F1);
f2=log(F2);
x1=log(X1);
x2=log(X2);
exp_vf = VF;

zap=log(Z); 
qqp=log(QQ); 
rp=log(R);
ptilp=log(PTIL); 
rkp=log(RK); 
mutilp=log(MUTIL); 
mcp=log(MC); 
sp=log(S); 
htilp=log(HTIL);
hdp = log(HD);
stilp=log(STIL); 
kp=log(K); 
hp=log(H); 
ivp=log(IV); 
wp=log(W); 
wtilp=log(WTIL); 
outputp=log(OUTPUT); 
gp=log(G); 
mhp=log(MH); 
cp=log(C); 
lap=log(LA); 
paip=log(PAI); 
up=log(U);
cbackp=log(C);
ivbackp=log(IV); 
rbackp=log(R); 
paibackp=log(PAI); 
outputbackp=log(OUTPUT); 
wbackp=log(W);
f1p=log(F1);
f2p=log(F2);
x1p=log(X1);
x2p=log(X2);
exp_vfp = VF;

%CEE_DEEP.M
%Deep Structural parameters of  the economy  described in ``Optimal Operational Monetary Policy in the Christiano-Eichenbaum-Evans Model of the U.S. Business Cycle,'' by S.\ Schmitt-Grohe and M.\ Uribe (2004). 
function [BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, SIGMAM, KAPA, CHI, CHITIL, GAMA1, GAMA2, PSSI, PHI0, PHI1, PAISTAR, RHOZ, RHOG, STD_EPSZ, STD_EPSG, Z, G, sig] = cee_deep
%This program produces the deep structural parameters and steady-state values of endogenous variables of a standard rbc model augmented to allow for price and wage  stickiness  
%(c) Stephanie Schmitt-Grohe and Martin Uribe, January 31, 2004. 


RHOZ = 0.979; %Serial correlation of tech. shock (King and Rebelo, Handbook of Macro, p. 953)

STD_EPSZ = 0.0072; %Std. dev.  of innovation to tech. shock (King and Rebelo, Handbook of Macro, p. 953)

RHOG = 0.96; %Serial correlation of gov't consumption (Christiano and Eichenbaum, AER Eq. (10) and table 1 

STD_EPSG = 0.02; %Std. Dev. of  innovation to gov't consumption (Christiano and Eichenbaum, AER Eq. (10) and table 1

HD=1; %Steady State labor demand

BETTA = 1.03^(-1/4); %Discount factor;

THETA = 0.36; %Capital share

DELTA = 0.025; %Depreciation rate

NU = 1; %Fraction of wage bill subject to working capital constraint

PAISTAR = 1.042^(1/4); %Inflation target

PAI = PAISTAR; %Inflation

ETATIL = 21; %Labor elasticity of subst

ETA = 6; %goods elasticity of substitution

ALFA =  0.6; %Degree of price stickiness

ALFATIL = 0.64; %degree of wage stickiness

B = 0.65; %degree of habit formation

SIGMAM = 10.62; %money demand parameter

U = 1; %SS capital utilization

GAMA2_O_GAMA1 = 0.01; %cost of capital utilization

KAPA = 2.48; %Capital adjustment cost

CHI = 1; %Degree of indexation

CHITIL = 1; %Degree of wage indexation

SG = 0.18; %Share of gov't expenditures in gdp (Christiano and Eichenbaum, AER


Z = 1; %SS productivity factor

QQ = 1; %Tobin's q

sig = 0.1797;%0.1; %parameter scaling the standard deviation of exogenous shocks

R = PAI / BETTA; %Nominal interest rate

PTIL = ((1-ALFA * PAI^((ETA-1)*(1-CHI))) / (1-ALFA)) ^ (1/(1-ETA)); %relative price of optimizing firms

WTIL_O_W = ((1-ALFATIL * PAI^((CHITIL-1)*(1-ETATIL))) / (1-ALFATIL)) ^ (1/(1-ETATIL)); 

RK = (1/BETTA - 1 + DELTA) * QQ; %rental rate of capital

MUTIL = (ETATIL / (ETATIL-1)) * ((1-ALFATIL * BETTA * PAI^((1-CHITIL)*(ETATIL-1))) / (1 - ALFATIL * BETTA * PAI^((1-CHITIL)*ETATIL))) * (1/WTIL_O_W); % inverse of wage markup

MC = PTIL * (ETA -1 ) * (1- ALFA * BETTA * PAI^(ETA*(1-CHI))) / (ETA * (1-ALFA*BETTA* PAI^((CHI-1)*(1-ETA)))); % marginal cost of producing a unit of final good

S = ((1-ALFA) * PTIL^(-ETA)) / (1-ALFA*PAI^((1-CHI) * ETA)); %distortion between output and total production

HTIL = WTIL_O_W^(-ETATIL); %hours supplied by those who can reset the wage this period      

STIL = (1-ALFATIL) * WTIL_O_W^(-ETATIL) / (1-ALFATIL * PAI^((1-CHITIL)*ETATIL));%distortion between hours supplied and hours used as factor input

GAMA1 = RK;%parameter of the cost of utilization function: a(u) = GAMA1 * (u-1) + GAMA2/2(u-1)^2

GAMA2 = GAMA2_O_GAMA1 * GAMA1; %%parameter of the cost of utilization function:

K = (RK/MC/THETA)^(1/(THETA-1)); %aggregate capital factor input

H = STIL; %number of hours supplied (> than those actually used at the firm level)

IV = DELTA * K; %investment

W = MC * (1-THETA) * K^THETA / (1 + NU*(1-1/R));%average real wage rate

WTIL = WTIL_O_W * W; %Wage charged by optimizing unions

PSSI = K^THETA - S*(RK*K + MC*(1-THETA)*K^THETA); %FIXED Cost

OUTPUT = (K^THETA - PSSI) / S; %aggregate demand = c + g +i 

C = OUTPUT*(1-SG)-IV;%Steady state consumption, one should subtract a(u)k, but this is zero

G = SG * OUTPUT; %Government spending

LA = (1-B*BETTA) / C / (1-B); %steady state marginal utility of income

PHI0 = LA * W / H / MUTIL; %preference parameter, weight on h_t^2 in the utility function

PHI1 =  0.5392564909446074; %preference parameter, weight on money in the utility function.

MH = ( LA *(1-BETTA/PAI)/PHI1)^(-1/SIGMAM); %Real balances held by households 

F1 = (LA * WTIL^ETATIL * HTIL / (1 - ALFATIL * BETTA * PAI^((1-CHITIL)*(ETATIL-1)))) / WTIL^ETATIL; 

F2 = (ETATIL-1) / ETATIL * WTIL * F1;

X2 = OUTPUT * PTIL^(-ETA) / (1-ALFA* BETTA * PAI^((CHI-1)*(1-ETA)));

X1 = (ETA-1) / ETA * X2;

UF = log(C*(1-B)) - PHI0/2 * H^2 + PHI1 * MH^(1-SIGMAM)/(1-SIGMAM);

VF = UF/(1-BETTA);

EXP_VF = exp(VF);

za=log(Z); 
qq=log(QQ); 
r=log(R);
ptil=log(PTIL); 
rk=log(RK); 
mutil=log(MUTIL); 
mc=log(MC); 
s=log(S); 
htil=log(HTIL);
hd = log(HD);
stil=log(STIL); 
k=log(K); 
h=log(H); 
iv=log(IV); 
w=log(W); 
wtil=log(WTIL); 
output=log(OUTPUT); 
g=log(G); 
mh=log(MH); 
c=log(C); 
la=log(LA); 
pai=log(PAI); 
u=log(U);
cback=log(C);
ivback=log(IV); 
rback=log(R); 
paiback=log(PAI); 
outputback=log(OUTPUT); 
wback=log(W);
f1=log(F1);
f2=log(F2);
x1=log(X1);
x2=log(X2);
exp_vf = log(EXP_VF);

zap=log(Z); 
qqp=log(QQ); 
rp=log(R);
ptilp=log(PTIL); 
rkp=log(RK); 
mutilp=log(MUTIL); 
mcp=log(MC); 
sp=log(S); 
htilp=log(HTIL);
hdp = log(HD);
stilp=log(STIL); 
kp=log(K); 
hp=log(H); 
ivp=log(IV); 
wp=log(W); 
wtilp=log(WTIL); 
outputp=log(OUTPUT); 
gp=log(G); 
mhp=log(MH); 
cp=log(C); 
lap=log(LA); 
paip=log(PAI); 
up=log(U);
cbackp=log(C);
ivbackp=log(IV); 
rbackp=log(R); 
paibackp=log(PAI); 
outputbackp=log(OUTPUT); 
wbackp=log(W);
f1p=log(F1);
f2p=log(F2);
x1p=log(X1);
x2p=log(X2);
exp_vfp = log(EXP_VF);


%CEE_MODEL.M
function [fx, fxp, fy, fyp, fypyp, fypy, fypxp, fypx, fyyp, fyy, fyxp, fyx, fxpyp, fxpy, fxpxp, fxpx, fxyp, fxy, fxxp, fxx, f]=cee_model(crule,brule,frule)

%This program computes the function f containing the equilibrium condiitons of the economy  described in ``Optimal Operational Monetary Policy in the Christiano-Eichenbaum-Evans Model of the U.S. Business Cycle,'' by S.\ Schmitt-Grohe and M.\ Uribe (2004).
%The function f defines  the DSGE model (a p denotes next-period variables) : 
%  E_t f(yp,y,xp,x) =0. 
%
%Inputs: (crule,brule,frule) defining the type of interest-rate rule (contemporaneous-, backward-, or forward-looking)
%
%Output: analytical expressions for f and its first and 2nd derivatives
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%
%Date February 2004


approx = 2; %Order of approximation desired

%Define parameters
syms BETTA THETA DELTA NU  ETATIL ETA ALFA ALFATIL B SIGMAM KAPA CHI CHITIL GAMA1 GAMA2 PSSI PHI0 PHI1 ALFA_PAI ALFA_Y ALFA_R PAISTAR RSTAR RHOZ RHOG OUTPUT G 

%Define variables 
%Notation: Changes w.r.t. the paper are iv for i and output for y. Also f1 and f2 here correspond to what in the paper are f1/wtil^ETATIL f2/wtil^ETATIL

syms za qq r ptil rk mutil mc s htil hd stil k h iv w wtil output g mh c la pai u f1 f2 x1 x2 exp_vf cback ivback rback paiback outputback wback 

syms zap qqp rp ptilp rkp mutilp mcp sp htilp hdp stilp kp hp ivp wp wtilp outputp gp cp lap paip up f1p f2p x1p  x2p exp_vfp cbackp ivbackp rbackp paibackp outputbackp wbackp  

%Notes
% xback = x_{t-1} and xbackp = x_t;
%change of notation w.r.t. paper: s = s_{t-1} and sp=s_t
%exp_vf = exp(value function). We introduce the exp of the v.f. rather than the level, because we want eventually to work with the level and this program converts all variables into logs below. 

%Production Function
pf = za* (u*k)^THETA * hd^(1-THETA) - PSSI;
pfp = zap * (up*kp)^THETA * hdp^(1-THETA) - PSSI;

%Capital utilization cost
a = GAMA1 * (u-1) + GAMA2/2 * (u-1)^2;
ap = GAMA1 * (up-1) + GAMA2/2 * (up-1)^2;

%Investment
iv  = output  - c  - g -  a*k;
ivp = outputp - cp - gp - ap*kp;

%Capital adjustment costs
qiv = iv * (1 - KAPA /2 * (iv/ivback -1)^2);

d1qiv = qiv/iv -KAPA*(iv/ivback-1)*iv/ivback;
d2qivp = (ivp/iv)^2*KAPA*(ivp/iv -1);

rk= diff(a,'u');
rkp = diff(ap,'up');

mh = ((la -  BETTA * lap / paip) / PHI1)^(-1/SIGMAM);

mc  = rk / (za  * THETA * (u *k /hd )^(THETA-1));
mcp = rkp /(zap * THETA * (up*kp/hdp)^(THETA-1));

w  = mc  * za  * (1-THETA) * (u  * k  / hd )^(THETA)   / (1 + NU * (r -1)/r );
wp = mcp * zap * (1-THETA) * (up * kp / hdp)^(THETA)   / (1 + NU * (rp-1)/rp );

ptil  = ((1 - ALFA * pai ^(ETA-1) * paiback ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));
ptilp = ((1 - ALFA * paip^(ETA-1) * paibackp^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

x1  = (ETA -1 ) / ETA * x2;
x1p = (ETA -1 ) / ETA * x2p;

wtil  = ((w ^(1-ETATIL) - ALFATIL * wback ^(1-ETATIL) * (paiback^CHITIL / pai) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));
wtilp = ((wp^(1-ETATIL) - ALFATIL * wbackp ^(1-ETATIL) * (paibackp^CHITIL / paip) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

htil = (wtil / w)^(-ETATIL) * hd;
htilp = (wtilp / wp)^(-ETATIL) * hdp;

h = stilp * hd;

mutil = la * w / (PHI0 * h);

f1 = - ETATIL * f2 / ((1 - ETATIL) * wtil);
f1p = - ETATIL * f2p / ((1 - ETATIL) * wtilp);

e1 = -kp + (1-DELTA) * k + qiv;

e2 = -la + 1 / (c-B*cback) - BETTA * B / (cp-B*c);

e4 = -la * qq + BETTA * lap * (rkp * up - ap + qqp * (1-DELTA));

e5 = -la + la * qq * d1qiv + BETTA * lap * qqp * d2qivp;  

e8 = -f1 + la  * htil + ALFATIL * BETTA * (paip / pai^CHITIL)^(ETATIL-1) * f1p* (wtilp/wtil)^ETATIL; 

e9 = -f2  + la / mutil * w^(1+ETATIL) * hd / wtil^ETATIL + ALFATIL * BETTA * (paip / pai^CHITIL)^ETATIL * f2p* (wtilp/wtil)^ETATIL;

%Interest-rate rule
e11 = -log(r/RSTAR) + ALFA_R * log(rback/RSTAR) + crule * (ALFA_PAI * log(pai/PAISTAR) + ALFA_Y * log(output/OUTPUT)) + brule * (ALFA_PAI * log(paiback/PAISTAR) + ALFA_Y * log(outputback/OUTPUT)) + frule * (ALFA_PAI * log(paip/PAISTAR) + ALFA_Y * log(outputp/OUTPUT));  

e12 = -la + BETTA * r * lap / paip;

e16 = -x1 + output * mc * ptil^(-ETA-1) + ALFA * BETTA * lap / la * (ptil / ptilp)^(-ETA-1) * (pai^CHI / paip)^(-ETA) * x1p;

e17 = -x2 + output * ptil^(-ETA) + ALFA * BETTA * lap / la * (pai^CHI / paip )^(1-ETA) * (ptil / ptilp)^(-ETA) * x2p;

e20 = -pf + (c+iv+g+a*k) * sp;

e21 = -sp + (1 - ALFA) * ptil^(-ETA) + ALFA * (pai/ paiback^CHI)^ETA * s;

e24 = -stilp + (1-ALFATIL) * (wtil/w)^(-ETATIL) + ALFATIL * (wback / w) ^(-ETATIL) * (pai / paiback^CHITIL)^ETATIL * stil;

%Value function
uf = log(c - B * cback) - PHI0 /2 * h^2 + PHI1 * mh^(1-SIGMAM)/ (1-SIGMAM);

e26 = -log(exp_vf) + uf + BETTA * log(exp_vfp);

%Government Purchases shock
e27 = log(gp/G) - RHOG * log(g/G);

%Evolution of Technolgy shock
e28 = log(zap) - RHOZ * log(za);

%Definition of past inflation
e29 = rbackp - r;

%Definition of past inflation
e30 = paibackp - pai;

%Definition of past output
e31 = outputbackp - output;

e32 = ivbackp - iv;

e33 = cbackp -c;

e34 = -wbackp + w;

%Create function f
f = [e1; e2; e4; e5; e8; e9;  e11; e12;  e16; e17; e20; e21; e24; e26; e27; e28; e29; e30; e31; e32; e33; e34];

% Define the vector of controls in period t, y, and t+1, yp, and the vector of states in period t, x, and t+1, xp

x = [cback ivback rback paiback outputback  s stil wback k g za];

xp = [cbackp ivbackp rbackp paibackp outputbackp  sp stilp wbackp kp gp zap];

y = [exp_vf qq r hd output c la pai u f2 x2];

yp= [exp_vfp qqp rp hdp outputp cp lap paip up f2p x2p];

%Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));

%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,x,y,xp,yp,approx);