// BAXTER INCOMPLETE MARKETS CASE, SIGMA=2, TEMP G SHOCK
var y, c, k, h, inve, p, b, g, tau, yStar, cStar, kStar, hStar, inveStar, pStar, bStar, taxRev, nx, nxStar;
varexo e_g, e_tau;

parameters  A, AStar, Gss, alpha, beta, delta, eta, gStar, gamma, pay, rho_g, rho_t, sigma, tauSS, tauStar, theta;
A = 1;
AStar = 1;
Gss = 0.15491;
alpha = 0.58;
beta  = 0.984;
delta = 0.025;
eta = 1; 
gStar = 0.15491;
gamma = 1.004;
pay = 0.83235;
rho_g = 0.90;
rho_t = 0.90;
sigma = 2;
tauSS  = 0.3;
tauStar = 0.3;
theta = 0.282;

model;
//c(+1)/c=cStar(+1)/cStar;

p= theta*(c^((theta*(1-sigma))-1))*((1-h)^((1-theta)*(1-sigma)));

pStar= theta*(cStar^((theta*(1-sigma))-1))*((1-hStar)^((1-theta)*(1-sigma)));

p(+1)/p=pStar(+1)/pStar;

1/(beta*c*(1-(inve/k(-1)-delta)/(eta*delta)))=
(1/(c(+1)*(1-(inve(+1)/k-delta)/(eta*delta))))*(1-delta+((inve(+1)/k)*(inve(+1)/k-delta)/(eta*delta))
-(1/2*eta*delta)*((inve(+1)/k-delta)^2))+(1/gamma)*(1-tau)*(1-alpha)*(y(+1)/k)*(1/c(+1));

1/(beta*cStar*(1-(inveStar/kStar(-1)-delta)/(eta*delta)))=
(1/(cStar(+1)*(1-(inveStar(+1)/kStar-delta)/(eta*delta))))*(1-delta+((inveStar(+1)/kStar)*(inveStar(+1)/kStar-delta)/(eta*delta))
-(1/2*eta*delta)*((inveStar(+1)/kStar-delta)^2))+(1/gamma)*(1-tauStar)*(1-alpha)*(yStar(+1)/kStar)*(1/cStar(+1));

gamma*k=(1-delta)*k(-1)+(inve/k(-1)-(1/2*eta*delta)*((inve/k(-1)-delta)^2))*k(-1);

gamma*kStar=(1-delta)*kStar(-1)+(inveStar/kStar(-1)-(1/2*eta*delta)*((inveStar/kStar(-1)-delta)^2))*kStar(-1);

(1/theta-1)*c/(1-h)=(1-tau)*alpha*y/h;

(1/theta-1)*cStar/(1-hStar)=(1-tauStar)*alpha*yStar/hStar;

y = A*(k(-1)^(1-alpha))*(h^alpha);

yStar = AStar*(kStar(-1)^(1-alpha))*(hStar^alpha);

pay*(y-c-inve-g)+(1-pay)*(yStar-cStar-inveStar-gStar)=0;

pay*b+(1-pay)*bStar=0;

bStar(-1)+yStar-cStar-inveStar-gStar-beta*(p(+1)/p)*bStar=0;

g = (1-rho_g)*Gss+rho_g*g(-1) + e_g;

tau = (1-rho_t)*tauSS+rho_t*tau(-1) + e_tau;

taxRev = tau*y;

nx = y-c-g-inve;

nxStar = yStar-cStar-gStar-inveStar;

//b(-1)+y-c-inve-g-beta*(p(+1)/p)*b=0;
//r=(gamma/beta)*(p(-1)/p)-1;

end;

initval;
y = 0.77453;
c = 0.49402;
h = 0.2;
k = 5.0239;
g = 0.15491;
tau = 0.3;
taxRev = 0.2323;
p = 0.81742;
b = 0;
nx = 0;
// r = 0.0203;
inve = 0.1256;
yStar = 0.77453;
cStar = 0.49402;
hStar = 0.2;
kStar = 5.0239;
inveStar = 0.1256;
bStar = 0;
pStar = 0.81742;
nxStar = 0;
e_g = 0;
e_tau = 0;
end;

check;

steady;

shocks;
var e_g; stderr 0.025946;
var e_tau; stderr 0;
//stderr 0.3892;
//sqrt(1.2973*0.3=0.3892)=0.06239
end;


stoch_simul(periods=1000, irf=40);
//order=1
