
%% Preliminaries
clear
addpath C:\Dynare\5.0.0\matlab

%% Calibrated parameters
n=2;

% Industry weights
Phi_1 = 0.5;
Phi_2 = 0.5;

% Consumption share: direct elasticity of the consumer price index Pc to changes in the price of the various goods
beta_1=0.5; % Sector 1
beta_2=0.5; % Sector 2

% Input-Ouput Matrix: direct elasticity of sector i's cost to a change in sector j's price
Omega_11=0.1;
Omega_12=0.3;
Omega_21=0.5;
Omega_22=0.2;

% Price adjusment parameters
delta_1=0.8; % Sector 1
delta_2=0.2; % Sector 2

% rho (dicount factor)
rho = 0.6;
par.rho = rho;

% Autoregression technology
rho_1 = 0.9;
par.rho_1 = rho_1;
rho_2 = 0.9;
par.rho_2 = rho_2;

% gamma
gamma=0.8;
par.gamma=gamma;

% phi
phi = 0.6;
par.phi=phi;

% zeta
zeta = 0.8;
par.zeta=zeta;

%% Build Parameters

% Consumption shared
beta = [beta_1 beta_2]'; par.beta = beta;
betaT = beta';
par.betaT_1 = betaT(1);
par.betaT_2 = betaT(2);

% Input-Ouput Matrix
Omega = [Omega_11,Omega_12;Omega_21,Omega_22]; par.Omega = Omega;

% Leontief (from Xi equation):
D1 = (eye(n)-Omega); par.D1 = D1;

% Price adjusment parameters
Delta_1 = (delta_1*(1-rho*(1-delta_1)))/(1-rho*delta_1*(1-delta_1));
Delta_2 = (delta_2*(1-rho*(1-delta_2)))/(1-rho*delta_2*(1-delta_2));
Delta = diag([Delta_1, Delta_2]); par.Delta = Delta;

% Domar weights (mnatural rate equation)
lambdaT=beta'*inv(eye(n)-Omega);
par.lambdaT_1 = lambdaT(1);
par.lambdaT_2 = lambdaT(2);

% labor share
alpha = (eye(n)-Omega)*[1;1]; 
par.alpha_1 = alpha(1);
par.alpha_2 = alpha(2);

% Aggregate Phillips Curve: coefficient on output
B = (Delta*inv(eye(n)-Omega*Delta)*alpha)/(1-beta'*Delta*inv(eye(n)-Omega*Delta)*alpha); par.B = B;

% Aggregate Phillips Curve: coefficient on inflation expectations
theta = (Delta*inv(eye(n)-Omega*Delta)-B*(lambdaT-beta'*Delta*inv(eye(n)-Omega*Delta)))*(eye(n)-Omega);
par.theta_11 = theta(1,1);
par.theta_12 = theta(1,2);
par.theta_21 = theta(2,1);
par.theta_22 = theta(2,2);

% From Sectoral Phillips Curve:
D2 = rho*(eye(n) - theta);
par.D2_11 = D2(1,1);
par.D2_12 = D2(1,2);
par.D2_21 = D2(2,1);
par.D2_22 = D2(2,2);

% From Sectoral Phillips Curve:
D3 = B*(gamma-phi);
par.D3_1 = D3(1);
par.D3_2 = D3(2);

% Leontief inverse: Inv(I_{n} - Omega)
D4 = inv(eye(n)-Omega);
par.D4_11 = D4(1,1);
par.D4_12 = D4(1,2);
par.D4_21 = D4(2,1);
par.D4_22 = D4(2,2);

% Industry weights
Phi = [Phi_1, Phi_2]'; par.Phi = Phi;

% Inflation aggregator
delta_phi = Phi'*Delta*inv(eye(n)-Omega*Delta);

% Inflation aggregator
delta_beta = beta'*Delta*inv(eye(n)-Omega*Delta);

% Aggregate slope of Phillips Curve
kappa_phi = (delta_phi*alpha)/(1-delta_beta*alpha); par.kappa_phi = kappa_phi;

% Inflation aggregator
D5 = rho*Phi'*(eye(n) - theta);
par.D5_1 = D5(1);
par.D5_2 = D5(2);

% Residual u of Phillips Curve
D6 = (delta_phi*alpha*((lambdaT-delta_beta)/(1-delta_beta*alpha))-delta_phi)*(eye(n) - Omega);
par.D6_1 = D6(1);
par.D6_2 = D6(2);

% For productivity shock
D7 = theta*inv(eye(n)-Omega);
par.D7_11 = D7(1,1);
par.D7_12 = D7(1,2);
par.D7_21 = D7(2,1);
par.D7_22 = D7(2,2);

% Law of motion of sectoral inflation rates
D8  = inv(eye(n)-Delta*Omega)*Delta*alpha;
par.D8_1 = D8(1);
par.D8_2 = D8(2);

D9 = (eye(n)-Omega);
par.D9_11 = D9(1,1);
par.D9_12 = D9(1,2);
par.D9_21 = D9(2,1);
par.D9_22 = D9(2,2);

D10 = inv(Delta)*(eye(n)-Delta*Omega);
par.D10_11 = D10(1,1);
par.D10_12 = D10(1,2);
par.D10_21 = D10(2,1);
par.D10_22 = D10(2,2);

D11 = inv(Delta)*rho*(eye(n)-Delta);
par.D11_11 = D11(1,1);
par.D11_12 = D11(1,2);
par.D11_21 = D11(2,1);
par.D11_22 = D11(2,2);

%% Run Dynare
par.quiet = 0; % Set par.quiet = 0 for IRFs and diagnostics
save parametervalues par % Saves par to read in Dynare
dynare model_a  % Runs model in Dynare                                       
