
clear all
close all

%% Calibrating the model

disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALIBRATING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')

Datainit = [1 1 2 3; 2 2 4 5; 3 3 0 0; 1 1 0 0; 4 4 0 0];

matsize = size(Datainit,2) - 2;
tii = size(Datainit,1) - 2;
VA = size(Datainit,1) - 1;
tio = size(Datainit,1);
ti = size(Datainit,2) - 1;
tup = size(Datainit,2);

invflows = [1 1; 1 1];
empcompensation = [0.5 0.5];

% Vector of consumer expenditure shares

psi = (Datainit(1:matsize,tup) - Datainit(1:matsize,ti))./sum(Datainit(1:matsize,tup) -...
    Datainit(1:matsize,ti));

% Vector of value added shares of firms

eta = Datainit(VA,1:matsize)./Datainit(tio,1:matsize);

% Employment share of firms

alpha = empcompensation./Datainit(VA,1:matsize);

% IO matrix

Omega = zeros(matsize,matsize);
for i = 1:matsize
    Omega(i,:) = Datainit(i,1:matsize)./Datainit(tii,1:matsize);
end
Omega(isnan(Omega)) = 0;

% Investment IO matrix

invIO = [0.5 0.5; 0.5 0.5];

% Share of sales to firm j out of total sales of a firm

Delta = zeros(matsize,matsize);
for i = 1:matsize
    Delta(:,i) = Datainit(1:matsize,i)./(Datainit(1:matsize,tup) +...
        sum(invflows,2));
end

% Share of investment good sales to firm j out of total sales of a firm

iota = zeros(matsize,matsize);
for i = 1:matsize
    iota(:,i) = invflows(1:matsize,i)./(Datainit(1:matsize,tup) +...
        sum(invflows,2));
end

% Share of final consumption sales out of the total sales of a firm

phi = (Datainit(1:matsize,tup) -...
    Datainit(1:matsize,ti))./(Datainit(1:matsize,tup) + sum(invflows,2));

% Value added share of each firm

lambda = Datainit(VA,1:matsize)./sum(Datainit(VA,1:matsize));

% Choosing parameters

demshock = ones(matsize,1);
demindex = 1; % Index of shocked sector
demshock(demindex) = 1.05; % Set demand shock

for i = 1:matsize
    if not(ismember(i,demindex))
        demshock(i) = (1-psi(demindex)'*(demshock(demindex)))/...
            (1-sum(psi(demindex)));
        % Normalizing shocks
    end
end

save("params.mat","alpha","eta","psi","Omega","invIO","lambda",...
    "Delta","iota","phi","demshock")

disp('Calibrated')

%% Solving model for effect of shock

disp('%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLVING MODEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')

% Calculating new steady state

% fun = @(variablesaux)eqsystem_fun_inv(variablesaux,alpha,psi,eta,invIO,iota,Omega,Delta,phi,lambda,demshock,matsize);
% 
% x0 = ones(2*matsize + 2,1);
% 
% solution = fsolve(fun,x0);
% 
% p_fut = solution(1:matsize); % Price hat
% output_fut = solution((matsize+1):(2*matsize)); % Output hat
% invprices_fut = exp(invIO'*log(p_fut)); % Invprice hat
% realrate_fut = invprices_fut; % Future real rates
% GDP_fut = solution(2*matsize+1); % GDP hat
% wage_fut = solution(2*matsize+2); % Wage hat
% k_fut = (p_fut.*output_fut)./realrate_fut; % Capital hat
% investment_fut = k_fut; % Investment hat
% 
% w = [ones(1,201) wage_fut]; 
% GDP = [ones(1,201) GDP_fut]; 
% p_1 = [ones(1,201) p_fut(1)];
% p_2 = [ones(1,201) p_fut(2)];
% r_1 = [ones(1,201) realrate_fut(1)];
% r_2 = [ones(1,201) realrate_fut(2)];
% y_1 = [ones(1,201) output_fut(1)];
% y_2 = [ones(1,201) output_fut(2)];
% I_1 = [ones(1,201) investment_fut(1)];
% I_2 = [ones(1,201) investment_fut(2)];
% k_1 = [ones(1,201) k_fut(1)];
% k_2 = [ones(1,201) k_fut(2)];
% 
% save("initvalues.mat","w","GDP","p_1","p_2","r_1","r_2","y_1","y_2","I_1",...
%     "I_2","k_1","k_2")

dynare dynprodnetmodel.mod