//Bansal-Yaron Long-Run Risks model

//original model
//x 		= rho*x(-1) + phi_e*vol(-1)*shock_x
//g 		= mu_c + x(-1) + vol(-1)*shock_g
//vol 	= (mu_s^2 + v*(((vol(-1))^2)-mu_s^2) + sigma_v*shock_s)^0.5;
//rf		= delta + (1-theta)*q + (mu_c+x(-1))/eis -0.5*(lambda_v^2*sigma_v^2) -((1-theta)*A2*(1-k1w*v)+0.5*(lambda_e^2+gamma^2))*vol(-1)^2
//r		= k0 + mu_d + (k1-1)*B0 +k1*B2*mu_s^2*(1-v) +B1*(k1*rho-1)*x(-1) + B2*(k1*v-1)*vol(-1) + k1*B1*phi_e*sqrt(vol(-1))*shock_x + k1*B2*sigma_v*shock_s


//variables
var r rf x g d vol;

//shocks
varexo shock_g shock_x shock_s shock_d;

//parameters
parameters delta gamma eis mu_c rho phi_e mu_s v sigma_v mu_d phi phi_d epsilon k0 k1 k0w k1w;

//initial parameter values
//fixed
k0 = -0.606855666;
k1 = 0.801052375;
k0w= 0.616217383; 
k1w= 0.182682555;

//to be estimated
delta = 0.03;
gamma = 2.5;
eis = 1.5;
mu_c = 0.01;
rho = 0.9;
phi_e = 0.2;
mu_s = 0.01;
v = 0.9;
sigma_v = 0.1;
mu_d =0.01;
phi = 5;
phi_d =0.2;
epsilon=0.5;

//model equations

model(linear);
x 		= rho*x(-1) + phi_e*vol(-1)*shock_x;
g 		= mu_c + x(-1) + vol(-1)*shock_g;
vol 	= (mu_s^2 + v*(((vol(-1))^2)-mu_s^2) + sigma_v*shock_s)^0.5;
d		= mu_d + phi*x(-1) + phi_d*vol(-1)*(epsilon*shock_g+(((1-epsilon^2)^(0.5))*shock_d));
rf		= delta + (1-((1-gamma)/(1-1/eis)))*(-delta + (1-1/eis)*mu_c + k0w + ((-delta+(1-1/eis)*mu_c+k0w+k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*mu_s^2*(1-v)+0.5*((1-gamma)/(1-1/eis))*(k1w^2)*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))^2*sigma_v^2)/(1-k1w))*(k1w-1) + k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*mu_s^2*(1-v)) + (mu_c+x(-1))/eis -0.5*(((1-((1-gamma)/(1-1/eis)))*k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v))))^2*sigma_v^2) - ((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*(phi_e^2)+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v) +0.5*(((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e)^2+gamma^2))*vol(-1)^2;
r		= k0 + mu_d + (k1-1)*(((v-1)*k1*(((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v)+0.5*gamma^2-gamma*phi_d*epsilon+0.5*phi_d +0.5*((k1*((phi-1/gamma)/(1-rho*k1))*phi_e)-((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e))^2)/(1-v*k1))*mu_s^2-0.5*((((1-gamma)/(1-1/eis))-1)*k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))+k1*(((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v)+0.5*gamma^2-gamma*phi_d*epsilon+0.5*phi_d +0.5*((k1*((phi-1/gamma)/(1-rho*k1))*phi_e)-((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e))^2)/(1-v*k1)))^2*sigma_v^2-(-delta + (((1-gamma)/(1-1/eis))-1)*(-delta + (1-1/eis)*mu_c + k0w + ((-delta+(1-1/eis)*mu_c+k0w+k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*mu_s^2*(1-v)+0.5*((1-gamma)/(1-1/eis))*(k1w^2)*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))^2*sigma_v^2)/(1-k1w))*(k1w-1) + k1w*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*mu_s^2*(1-v)) - mu_c/eis + k0 + mu_d))/(k1-1)) +k1*(((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v)+0.5*gamma^2-gamma*phi_d*epsilon+0.5*phi_d +0.5*((k1*((phi-1/gamma)/(1-rho*k1))*phi_e)-((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e))^2)/(1-v*k1))*mu_s^2*(1-v) +((phi-1/gamma)/(1-rho*k1))*(k1*rho-1)*x(-1) + (((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v)+0.5*gamma^2-gamma*phi_d*epsilon+0.5*phi_d +0.5*((k1*((phi-1/gamma)/(1-rho*k1))*phi_e)-((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e))^2)/(1-v*k1))*(k1*v-1)*vol(-1) + k1*((phi-1/gamma)/(1-rho*k1))*phi_e*sqrt(vol(-1))*shock_x + k1*(((1-((1-gamma)/(1-1/eis)))*(((1-gamma)/(1-1/eis))*((k1w^2)*((1-1/eis)/(1-k1w*rho))^2*phi_e^2+(1-1/eis)^2)/(2*(1-k1w*v)))*(1-k1w*v)+0.5*gamma^2-gamma*phi_d*epsilon+0.5*phi_d +0.5*((k1*((phi-1/gamma)/(1-rho*k1))*phi_e)-((1-((1-gamma)/(1-1/eis)))*k1w*((1-1/eis)/(1-k1w*rho))*phi_e))^2)/(1-v*k1))*sigma_v*shock_s;
end;

shocks;
var shock_x; stderr .2;
var shock_g; stderr .3;
var shock_s; stderr .1;
var shock_d; stderr .1;
end;



estimated_params;
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
// inverse-gamma for variances, beta for [0,1], normal for the rest.

delta	, 0.03	, 0.0001, 0.1	, BETA_PDF		, 0.03 , 0.1;
gamma	, 2.5	, 0.5	, 50	, NORMAL_PDF	, 2.5 , 10;
eis		, 1.5	, 0.0001, 5		, NORMAL_PDF	, 1.5 , 1;

mu_c	, 0.01	, 0.0001, 0.1000, NORMAL_PDF	, 0.01 , 0.010;
mu_d	, 0.01	, 0.0001, 0.0100, NORMAL_PDF	, 0.01 , 0.011;
mu_s	, 0.001	, .00001, 0.0100, NORMAL_PDF	, 0.01 , 0.011;
rho		, 0.99	, 0.9	, 0.9999, BETA_PDF		, 0.99 , 0.010;

phi		, 5		, 0.0001, 10	, NORMAL_PDF	, 5   , 5;
phi_e	, 0.2	, 0.0001, 0.9999, NORMAL_PDF	, 0.2 , 1;
phi_d	, 0.2	, 0.0001, 0.9999, NORMAL_PDF	, 0.2 , 1;
v		, 0.9	, 0.0001, 0.9999, BETA_PDF		, 0.9 , 0.1;
sigma_v	, 0.005 , 0.005 , 0.9999, INV_GAMMA_PDF	, .005, .005;
epsilon	, 0.5	, 0.0001, 0.9999, NORMAL_PDF	, 0.5 , 0.5;
end;

varobs g d rf r;

//bayes
estimation(datafile=LLR,nograph,xls_sheet=ready,xls_range=g2:j121,first_obs=20,mh_replic=5000,mh_nblocks=2);

//mle
//estimation(datafile=LLR,xls_sheet=ready,xls_range=g2:j121,nograph,mode_compute=6,first_obs=10,mh_replic=0,mh_nblocks=0,mh_jscale=0.20,mh_drop=0.2);