//fully deterministic, set betat and ystar fixed
//same as sgu_stoch_beta, but make A exogenous but deterministic variable
//path of A fully known and expected, ystar and foreign beta still stochastic
//do transitional dynamics, find out path as A changes over time
//make foreign discount factor betat stochastic, follow AR(1) process
//utility specification as in SGU

close all;

//define 16 unknowns, 3 lg multipliers and 2 shocks
var c, l, k, i, b, ch, cf, bcl, %home
    cstar, bstar, cstarf, cstarh, %foreign
    pf, p, pstar, rstar, %clearing 
    lambda, eta, %2 home lg multipliers
    nu, %foreign lg multiplier
    /*ystar, betat, btilde, %stochastic shocks*/
    y, byratio; %for checking results

varexo A /*, epsilon_y, epsilon_b*/; %kill shocks to ystar and betat
//A still exogenous but path known and expected, not stochastic

//define parameters
parameters alpha, phi, delta, sigma, chi, kappa, gamma, rho, gammastar, beta, psi, omega, N,
           zeta, tau, %allow coefficient and exponent on labor to be different between utility and discounting
           mu, xi, ind, %indicator for labor in discounting
           mu_a, mu_y, mu_b, a_ss, y_ss, b_ss, sigma_a, sigma_y, sigma_b, ystar;
           %these stochastic shocks don't apply here, ignore

alpha = 1/3; %weight of capital in production
phi = 4.5; %4.5 or 8, capital adjustment cost, changing phi does not change steady state at all, and not irf by much
delta = 0.025; %0.025, 0.1 in SGU, capital depreciation rate
sigma = 2; %exponent in utility
chi = 1; %1, 0, does not apply here
kappa = 1; %1, 2.95, disutility of labor, does not apply here

rho = 1 - 1/0.9; %elasticity of substitution is 0.9 or 1.5
gamma = 0.7; %ratio of home goods in home consumption, decreasing, 0.35 - 0.66 by world bank data
gammastar = 0.95; %ratio of foreign good in ROW consumption, decreasing, 0.95 - 0.99 by world bank data
                  //0.78 and 0.95, 0.8 and 0.961
psi = 0.11; %0.11, internal discounting exponent from SGU
omega = 1.455; %1.455, internal discounting factor from SGU
tau = 1; %1 is same as in SGU discounting
zeta = omega; %equal to omega is similar to SGU discounting
mu = 1.057; %1 is same as in SGU discounting, >1.085 gives positive b for rho 0.9
            %graphs look better and more robust than ghh?
            %1.085-1.09 for rho 0.9, 1.141-1.145 for rho 1.5, initial l is 2 
xi = 1; %1 is same as in SGU discounting
ind = 1; %when including labor in beta(c,l), 0 when dropping it

//changed psi,omega,beta to make b positive
beta = 0.98; %foreign discount factor
N = 4.2; %ratio between ROW population and China's
mu_a = 0.95;
mu_y = 0.95;
mu_b = 0.95;
a_ss = 1; %SS productivity factor
y_ss = 1; %SS ratio between ROW GDP and China GDP, 1.835
b_ss = 1; %in steady discount factor state equals beta
sigma_a = 0.01;%
sigma_y = 0.01; %repress irfs to y* shocks, don't care about it
sigma_b = 0.01;
//above all useless here, ystar always fixed, equal to 1
ystar = y_ss;

//system of raw equilibrium conditions with expectation and lg multipliers
model;
//%home FOC's
(c-tau*l^zeta/zeta)^(-sigma) - lambda*p=0; %1
-lambda + eta*( 1 - phi*(i/k(-1)-delta) ) = 0; %2
-lambda*pf + bcl*rstar*lambda(+1)*pf(+1) = 0; %3
-eta + bcl*( lambda(+1)*A(+1)*alpha*k^(alpha-1)*l(+1)^(1-alpha) 
+ eta(+1)*((1-delta)-0.5*phi*(delta-i(+1)/k)*(delta+i(+1)/k)) ) = 0; %4
-(c-tau*l^zeta/zeta)^(-sigma)*tau*l^(zeta-1) + lambda*(1-alpha)*A*k(-1)^alpha*l^(-alpha) = 0; %5

%foreign FOC's, no more betat
1/cstar - nu*pstar = 0; %6
-nu + beta*rstar*nu(+1) = 0; %7

%home budget
p*c + i + pf*b = A*k(-1)^alpha*l^(1-alpha) + rstar(-1)*pf*b(-1); %8
k = (1-delta)*k(-1) + i - phi*0.5*(i/k(-1)-delta)^2*k(-1); %9
%foreign budget
pstar*cstar + bstar = ystar + rstar(-1)*bstar(-1); %10

%bonds and good market clearing
b = -N*bstar; %11
ch + N*cstarh + i = A*k(-1)^alpha*l^(1-alpha); %12

%Uzawa internal discounting like is SGU
bcl = mu*(xi+c-ind*l^omega/omega)^(-psi); %13

%static cost minimization solution
p = ( gamma + (1-gamma)*pf^(rho/(rho-1)) )^((rho-1)/rho); %14
ch = c*gamma*p^( 1/(1-rho) ); %15
cf = c*(1-gamma)*(pf/p)^( 1/(rho-1) ); %16
pstar = ( gammastar + (1-gammastar)*pf^(rho/(1-rho)) )^((rho-1)/rho); %17
cstarf = cstar*gammastar*pstar^( 1/(1-rho) ); %18
cstarh = cstar*(1-gammastar)*(pf*pstar)^( 1/(1-rho) ); %19

%kill endogenous foreign discounting factor
//betat = beta*btilde;
%remove evolution of 3 shocks, each logged to ensure positivity
//A deterministically exogenous, no more endogenous, remove its transition
//ln(A) - ln(a_ss) = mu_a*( ln(A(-1)) - ln(a_ss) ) + epsilon_a; %20
//ln(ystar) - ln(y_ss) = mu_y*( ln(ystar(-1)) - ln(y_ss) ) + epsilon_y; %21
%minus sign in front of shock to see case of foreigners becoming impatient
//ln(btilde) - ln(b_ss) = mu_b*( ln(btilde(-1)) - ln(b_ss) ) - epsilon_b;

%for checking results with data
y= A*k(-1)^alpha*l^(1-alpha);
byratio = b/y; %should be close to 0.5 in real world data

end;

//define initial low A of economy and initial values for computation
/*initval;
A=0.9;
k=17;
rstar=0;
b=0;
bstar= -b/N;
end;*/

initval;

A = 0.9; %initial low productivity
//ystar = y_ss;
//btilde = b_ss;
//b = beta*btilde;
rstar = 1/beta;
bcl = 1/rstar;
k = 17;
l = k / ( (1/bcl-1+delta)/(alpha*A) )^(1/(alpha-1));%determined by initial k
c = (bcl/mu)^(-1/psi) - xi + ind*l^omega/omega;
i = delta*k;
p = (1-alpha)*A*k^alpha*l^(-alpha-zeta+1)/tau;
pf = ( (p^(rho/(rho-1)) - gamma) / (1-gamma) )^( (rho-1)/rho );
b = ( p*c+i-A*k^alpha*l^(1-alpha) )/((rstar-1)*pf);
bstar = -b/N;
pstar = ( gammastar + (1-gammastar)*pf^(rho/(1-rho)) )^((rho-1)/rho);
ch = c*gamma*p^( 1/(1-rho) );
cf = c*(1-gamma)*(pf/p)^( 1/(rho-1) );
cstar = (ystar+(rstar-1)*bstar) / pstar;
cstarf = cstar*gammastar*pstar^( 1/(1-rho) );
cstarh = cstar*(1-gammastar)*(pf*pstar)^( 1/(1-rho) );
lambda = (c-tau*l^zeta/zeta)^(-sigma)/p;
eta = lambda / ( 1 - phi*(i/k-delta) );
nu = 1 / (pstar*cstar);
y= A*k^alpha*l^(1-alpha);
byratio = b/y;

end;
steady(solve_algo=3);

//define ending state of high A 
endval;

A = 1; %ending with high productivity
//ystar = y_ss;
//btilde = b_ss;
//b = beta*btilde;
rstar = 1/beta;
bcl = 1/rstar;
l=1.5;
c = (bcl/mu)^(-1/psi) - xi + ind*l^omega/omega;
k = ( (1/bcl-1+delta)/(alpha*A) )^(1/(alpha-1))*l;
i = delta*k;
p = (1-alpha)*A*k^alpha*l^(-alpha-zeta+1)/tau;
pf = ( (p^(rho/(rho-1)) - gamma) / (1-gamma) )^( (rho-1)/rho );
b = ( p*c+i-A*k^alpha*l^(1-alpha) )/((rstar-1)*pf);
bstar = -b/N;
pstar = ( gammastar + (1-gammastar)*pf^(rho/(1-rho)) )^((rho-1)/rho);
ch = c*gamma*p^( 1/(1-rho) );
cf = c*(1-gamma)*(pf/p)^( 1/(rho-1) );
cstar = (ystar+(rstar-1)*bstar) / pstar;
cstarf = cstar*gammastar*pstar^( 1/(1-rho) );
cstarh = cstar*(1-gammastar)*(pf*pstar)^( 1/(1-rho) );
lambda = (c-tau*l^zeta/zeta)^(-sigma)/p;
eta = lambda / ( 1 - phi*(i/k-delta) );
nu = 1 / (pstar*cstar);
y= A*k^alpha*l^(1-alpha);
byratio = b/y;

end;
steady(solve_algo=3);

/*
//specify path of A
shocks;
//set path for exogenous deterministic A
var A;
periods 1 2 3 4 5 6 7 8 9;
values 0.991 0.992 0.993 0.994 0.995 0.996 0.997 0.998 0.999; 

//kill nexogenous stochastic shocks
//var epsilon_y = sigma_y^2;
//var epsilon_b = sigma_b^2;
end;
*/

//transitional dynamics for a set path of A
simul(periods=300);
rplot c;
rplot cstar;
rplot l;
rplot i;
rplot k;
rplot pf;
rplot rstar;
rplot bcl;
rplot b;


//can't do under deterministic case impulse response graphs for schocks to betat and ystar
//stoch_simul(order=1,irf=50) c,l,k,i,b,cstar,bstar,pf,rstar,bcl,A,ystar,betat;