 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Create VGR time series of real
% Y, C, G, I, NEx, Depreciation rate, Capital stock, labour & Years Schooling
% for GDR 1950/1960/1970-1989
%
% Daniel Fehrle
% June 2023,
% update: 12.12.2023 GDR depreciation rates, federal archive
% update: 04.06.2025 GDR CPI growth rates and r
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
net=1;
%%%%%%%%%%%%%% good quantities and stocks 
% based on Heske deflators, incl. sub aggregates price deflators


vgr_current=readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'Exp, current prices','Range','B2:M21');
vgr_100=readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'Exp, 1950=100','Range','B2:J41');
pop=readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'PopEmp','Range','B2:C37');


depr_pre70=readmatrix("depr.csv");
mps_net_inv=(depr_pre70(2,:)+depr_pre70(3,:))';
mps_net_inv=mps_net_inv./mps_net_inv(end); %%%==Heske Bruttoinv. nominal
depr_pre70=depr_pre70(2,:)';
depr_pre70=depr_pre70(1:end)./depr_pre70(end);



%deflate subaggregates
vgr_current_use=vgr_current(:,[3 4 6 7 8 9]);
vgr_100_use=vgr_100(:,[3 4 6 7 8 9])./vgr_100(end,[3 4 6 7 8 9]);

vgr_1989_use=vgr_current_use(end,:).*vgr_100_use;

%deflate depr. 
Pi=vgr_current_use(:,3)./vgr_1989_use(21:end,3);
depr_1989=vgr_current(:,11)./Pi;
depr_1989=fillmissing(depr_1989,'spline');

Pi_pre70=mps_net_inv./(vgr_100(11:21,[6])./vgr_100(21,[6]));
depr_pre70=depr_pre70./Pi_pre70;

depr_1989=[depr_1989(1)*depr_pre70(1:end-1); depr_1989];

depr_1989=[NaN(10,1); depr_1989];
%%%%%%%%%%%aggregates
%calculate GDP
gdp_1989=sum(vgr_1989_use(:,1:end-1),2)-vgr_1989_use(:,end);

%VA
vaf=vgr_current(:,10)./vgr_current(:,12);
va_current=[NaN(20,1);vaf];
va_1989=gdp_1989.*(exp(fillmissing(log(va_current-1),'linear'))+1);

%NEX
Nex_1989=vgr_1989_use(:,5)-vgr_1989_use(:,6);
du_1989=gdp_1989-Nex_1989;

vgr_1989=[vgr_1989_use Nex_1989 gdp_1989 va_1989 du_1989];

vgr_1989=fillmissing(vgr_1989,'spline');

%capital stock
%see Source excel
if net==0
capital_1990=1697332.15222455; %Gross
else
capital_1990=1211261.99014326; %Net
end




capital_1989(40)=capital_1990-vgr_1989(end,3)+depr_1989(end);
capital_1989=capital_1989';

for i=1:29

    capital_1989(end-i)=capital_1989(end+1-i)-vgr_1989(end-i,3)+depr_1989(end-i);
        
end


rate=depr_1989./capital_1989;
rate=fillmissing(rate,'linear');

for i=30:39

    capital_1989(end-i)=(capital_1989(end+1-i)-vgr_1989(end-i,3))./(1-rate(end-i)); 
        
end



depr=-diff(capital_1989)+vgr_1989(1:end-1,3);
depr_1989=[depr; depr_1989(end)];

rate_p70=mean(depr_1989(21:end)./capital_1989(21:end),'omitnan');



vgr_1989=[vgr_1989 depr_1989 capital_1989];

%%%%%%%%%%%%%%%%%%%%%%%%%%% CPI

%%nominal from Destatis 1970 -- 1989
nom_cons=fillmissing(vgr_current_use(:,1),'spline');
nom_cons=nom_cons(:,1)./nom_cons(end,1); %%%%norm on last observation
real_cons=vgr_1989(21:end,1)./vgr_1989(end,1); %%%%norm on last observation
CPI=nom_cons./real_cons;
pi_c=CPI(2:end)./CPI(1:end-1); 

%%nominal from MPS SZS 1960 -- 1970

nom_cons=readmatrix("MPS_cons.xls","Range","B2:B12");
nom_cons=nom_cons(:,1)./nom_cons(end,1); %%%%norm on last observation
real_cons=vgr_1989(11:21,1)./vgr_1989(21,1); %%%%norm on last observation
CPI=nom_cons./real_cons;
pi_c=[NaN;CPI(2:end)./CPI(1:end-1); pi_c];

%%%%%%%%%%%%%%%%%%%%%%%%%% r

r=[ones(11,1)*1.03;ones(19,1)*1.0325];


%%%%money market

MM=[(1960:1989)',pi_c,r];

writematrix(MM,'mny_e.csv');



%%%%%%%%%%%%%%%%%%%%%%hours worked, employed

lf=fillmissing(readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'PopEmp','Range','C3:C42'),'spline');
hours=fillmissing(readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'Alt_real GDP, Hours','Range','B2:B41'),'linear');
hwrk=hours.*lf;

%%%%%%%%%%%%%%%%%%%%%%years schooling
tyr=readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'Human','Range','B2:C42');
tyr(21:end,:)=fillmissing(tyr(21:end,:),'spline');
tyr=fillmissing(tyr,'linear');
tyr_pc=tyr(1:end-1,:);



%per capita

pop=fillmissing(readmatrix("Source_NIPA_EAST_GER.xls","Sheet",'PopEmp','Range','B3:B42'),'spline');
tyr=tyr_pc.*pop; 
vgr_1989_pc=vgr_1989./pop;
hwrk_pc=hwrk./pop;
lf_pc=lf./pop;


%
names={'T','C','G','I','V','Ex','Im','Nex','GDP','VA','DU','Depr','Capital','Labor Force','H','YS25','YS15','pop'};
names_pc={'T','c','g','i','v','ex','im','nex','gdp','va','du','depr','capital','labor force','h','ys25','ys15','pop'};
data=array2table([(1950:1989)' vgr_1989 lf hwrk tyr pop],'VariableNames',names);
data_pc=array2table([(1950:1989)' vgr_1989_pc lf_pc hwrk_pc tyr_pc pop],'VariableNames',names_pc);

if net==0
writetable(data,'abs.csv','WriteRowNames',true);
writetable(data_pc,'pc.csv','WriteRowNames',true);
else
writetable(data,'abs_net.csv','WriteRowNames',true);
writetable(data_pc,'pc_net.csv','WriteRowNames',true);
end


%%%%%%%%%%%%%%per BIP for comparison

vgr_BIP=[vgr_1989(:,[1 2 3]) vgr_1989(:,4)+vgr_1989(:,7)]./vgr_1989(:,8);
cap_BIP=vgr_1989(:,end)./vgr_1989(:,8);
depr_cap=vgr_1989(:,end-1)./vgr_1989(:,end);
GDP_growth=[NaN; vgr_1989(2:end,8)./vgr_1989(1:end-1,8)];
pop_growth=[NaN; pop(2:end)./pop(1:end-1)];
names_BIP={'T','cy','gy','iy','nexy','yg', 'H', 'tyr', 'gpop'};
data_BIP=array2table([(1950:1989)' vgr_BIP GDP_growth hwrk_pc tyr_pc(1:end,2) pop_growth cap_BIP depr_cap]);

if net==0
writetable(data_BIP,'py.csv','WriteRowNames',false);
else
writetable(data_BIP,'py_net.csv','WriteRowNames',false);
end

