
addpath C:\dynare\6.1\matlab
clear all
close all

%% Calibrating the model

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

OmegaRaw = [5 5; 7 7]; % Input-output flow matrix
finaluse = [20; 20]; % Final use vector
Govtsec = [2;2]; % Government sector
invsales = [4;4]; % Investment sales
invflows = [2 2; 2 2];

disc = 0.98;

matsize = size(OmegaRaw,1);

% Employment share of firms

alpha = [0.5 0.5];

% Depreciation rates

py = finaluse + sum(OmegaRaw,2);
depr = zeros(2,1);
for i = 1:matsize
    depr(i) = ((1-disc)/disc)*(invsales(i)/...
        ((py(i) - sum(OmegaRaw(:,i),1))*(1-alpha(i))))/(1-(invsales(i)/...
        ((py(i) - sum(OmegaRaw(:,i),1))*(1-alpha(i)))));
end

% Vector of value added shares of firms

eta = 1 - sum(OmegaRaw,1)./(sum(OmegaRaw,2) + finaluse)';

% IO matrix

Omega = zeros(matsize,matsize);
for i = 1:matsize
    Omega(i,:) = OmegaRaw(i,:)./sum(OmegaRaw,1);
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) = OmegaRaw(:,i)./py;
end

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

iota = zeros(matsize,matsize);
for i = 1:matsize
    for j = 1:matsize
        iota(i,j) = invIO(i,j)*...
            (disc*depr(j)/(1 - disc*(1-depr(j))))*...
            (1-alpha(j))*eta(j)*(py(j)/py(i));
    end
end

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

phi = (finaluse - Govtsec - invsales)./py;

% Share of government demand out of the total sales of a sector

gamma = Govtsec./py;

% Vector of consumer expenditure shares

psi = (finaluse - invsales - Govtsec)./sum(finaluse - invsales - Govtsec);

% Labour share out of total of each sector

lambda = (alpha.*(py' - sum(OmegaRaw,1)))./sum((alpha.*(py' - sum(OmegaRaw,1))));

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

disp('Calibrated')

%% Solving model for effect of shock

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

dynare dynprodnetmodel.mod
