% Read in Value of capital stock and smooth to create a monthly series from 1980:M1 through 2023:M6

clear all;
small = 1.0e-10;
pinv_tol = 1.0e-05;
big = 1.0e+8;

% -- File Directories   
outdir = 'out/';
figdir = 'fig/';
matdir = 'mat/';
datadir = '../Data/WeatherDamages/';


% Read in Capital Stock Data -- BEA Data
% Annual data from 1979 through 2022
str_data = [datadir 'BEA_NetStock_Table.xlsx'];
bea_dat = readmatrix(str_data,'Range','L8:BC8')';
yr_bea = (1979:2022)';
plot(yr_bea,bea_dat)

% Interpolate to create a monthly series from 1980:M1 through 2023:M6
% Use the PWT data for 1970 through 2019
% Use the BEA data for 2020 through 2023
t_vec = (1979:2022);
t_vec = t_vec+0.5;
t_vec = t_vec';
ln_bea = log(bea_dat);

% Interpolate and extrapolate
calvec = (1980:1/12:2023+5/12)';
Y = ln_bea;
tx = t_vec(1:size(Y,1),1);
% Dates for interpolation
t_int = (1980:0.083333333333333:tx(end))';
% Interpolate
Y_int_bea = csapi(tx,Y,t_int);
Y_bea = Y;
% Extrapolate
n_ext = size(calvec,1)-size(t_int,1);
% average value last 12 months
tmp = dif(Y_int_bea,1);
m = mean(tmp(end-11:end));
% Extrapolate
Y_ext = Y_int_bea(end)*ones(n_ext,1);
Y_ext = Y_ext + m*(1:n_ext)';
Y_int_bea = [Y_int_bea;Y_ext];

plot(calvec,Y_int_bea,tx,Y_bea,'o');

% Save These data series
KStock_BEA = exp(Y_int_bea);
KStock_BEA = KStock_BEA/KStock_BEA(1);


save([matdir 'Kstock.mat'],'calvec','KStock_BEA');   

