
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) = 2; % 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 military shock

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

dynare dynprodnetmodel.mod