%--------------------------------------------------------------------------
% Collateral Shocks
% Yvan Becard and David Gauthier, June 2020
% Data for Estimation
%--------------------------------------------------------------------------

clear

%--------------------------------------------------------------------------
% Load Data
%--------------------------------------------------------------------------

% Path
Path = '.\';

% Time format
timeformat = 'mm/dd/yyyy';

% Time span
startyear = 1985; endyear = 2019;
timespan = timeseries(ones,datenum(datetime(startyear,01,01):calmonths(3):datetime(endyear,01,01))');


%--------------------------------------------------------------------------
% Macroeconomic Time Series
%--------------------------------------------------------------------------

% Gross Domestic Product, Billions of Dollars, Quarterly, Seasonally Adjusted
[gdpraw,~,all] = xlsread([Path,'GDP.xls']);
gdpraw         = timeseries(gdpraw,datenum(all(12:end,1),timeformat));

% Gross Domestic Product: Implicit Price Deflator, Index 2009=100, Quarterly, Seasonally Adjusted
[gdpdeflatorraw,~,all] = xlsread([Path,'GDPDEF.xls']);
gdpdeflatorraw         = timeseries(gdpdeflatorraw,datenum(all(12:end,1),timeformat));

% Personal Consumption Expenditures: Nondurable Goods, Billions of Dollars, Quarterly, Seasonally Adjusted
[consumptionnondurableraw,~,all] = xlsread([Path,'PCND.xls']);
consumptionnondurableraw         = timeseries(consumptionnondurableraw,datenum(all(12:end,1),timeformat));

% Personal Consumption Expenditures: Services, Billions of Dollars, Quarterly, Seasonally Adjusted
[consumptionservicesraw,~,all] = xlsread([Path,'PCESV.xls']);
consumptionservicesraw         = timeseries(consumptionservicesraw,datenum(all(12:end,1),timeformat));

% Personal Consumption Expenditures: Durable Goods, Billions of Dollars, Quarterly, Seasonally Adjusted
[consumptiondurableraw,~,all] = xlsread([Path,'PCDG.xls']);
consumptiondurableraw         = timeseries(consumptiondurableraw,datenum(all(12:end,1),timeformat));

% Gross Private Domestic Investment, Billions of Dollars, Quarterly, Seasonally Adjusted
[investmentraw,~,all] = xlsread([Path,'GPDI.xls']);
investmentraw         = timeseries(investmentraw,datenum(all(12:end,1),timeformat));

% Gross private domestic investment (implicit price deflator), Index 2012=100, Quarterly, Seasonally Adjusted
[investmentdeflatorraw,~,all] = xlsread([Path,'A006RD3Q086SBEA.xls']);
investmentdeflatorraw         = timeseries(investmentdeflatorraw,datenum(all(12:end,1),timeformat));

% Nonfarm Business Sector: Hours of All Persons, Index 2009=100, Quarterly, Seasonally Adjusted
[hoursraw,~,all] = xlsread([Path,'HOANBS.xls']);
hoursraw         = timeseries(hoursraw,datenum(all(12:end,1),timeformat));

% Effective Federal Funds Rate, Percent, Quarterly, Average, Not Seasonally Adjusted
[fedfundsrateraw,~,all] = xlsread([Path,'FEDFUNDS.xls']);
fedfundsrateraw         = timeseries(fedfundsrateraw,datenum(all(12:end,1),timeformat));

% Civilian Noninstitutional Population, Thousands of Persons, Quarterly, Average, Not Seasonally Adjusted
[populationraw,~,all] = xlsread([Path,'CNP16OV.xls']);
populationraw         = timeseries(populationraw,datenum(all(12:end,1),timeformat));


%--------------------------------------------------------------------------
% Financial Time Series
%--------------------------------------------------------------------------

% Households and nonprofit organizations; debt securities and loans; liability, Level, Billions of Dollars, Quarterly, Seasonally Adjusted
[HNcreditraw,~,all] = xlsread([Path,'CMDEBT.xls']);
HNcreditraw         = timeseries(HNcreditraw,datenum(all(12:end,1),timeformat));

% 30-Year Fixed Rate Mortgage Average in the United States, Percent, Quarterly, Not Seasonally Adjusted
[HNrateraw,~,all] = xlsread([Path,'MORTGAGE30US.xls']);
HNrateraw         = timeseries(HNrateraw,datenum(all(12:end,1),timeformat));

% Nonfinancial business; debt securities and loans; liability, Level, Billions of Dollars, Quarterly, Seasonally Adjusted
[NFcreditraw,~,all] = xlsread([Path,'TBSDODNS.xls']);
NFcreditraw         = timeseries(NFcreditraw,datenum(all(12:end,1),timeformat));

% Moody's Seasoned Baa Corporate Bond Yield, Percent, Quarterly, Not Seasonally Adjusted
[NFrateraw,~,all] = xlsread([Path,'BAA.xls']);
NFrateraw         = timeseries(NFrateraw,datenum(all(12:end,1),timeformat));

%--------------------------------------------------------------------------
% Oil Time Series
%--------------------------------------------------------------------------
% Global price of brent crude; US Dollars per barrel, Dollars, Quarterly, Not Seasonally Adjusted
[Oilpriceraw,~,all] = xlsread([Path,'POILBREUSDQ.xls']);
Oilpriceraw         = timeseries(Oilpriceraw,datenum(all(12:end,1),timeformat));

% U.S field production of crude oil (thousand barrels)
[Oilproductionraw,~,all] = xlsread([Path,'PET_CRD_CRPDN_ADC_MBBL_M.xls']);
Oilproductionraw = timeseries(Oilproductionraw,datenum(all(12:end,1),timeformat));

%--------------------------------------------------------------------------
% Harmonize Units
%--------------------------------------------------------------------------

gdpraw                   = gdpraw                   * 1e9;
consumptionnondurableraw = consumptionnondurableraw * 1e9;
consumptionservicesraw   = consumptionservicesraw   * 1e9;
consumptiondurableraw    = consumptiondurableraw    * 1e9;
investmentraw            = investmentraw            * 1e9;
populationraw            = populationraw            * 1e3;
HNcreditraw              = HNcreditraw              * 1e9;
NFcreditraw              = NFcreditraw              * 1e9;
Oilproductionraw         = Oilproductionraw         * 1e3; 

%--------------------------------------------------------------------------
% Smooth Population Series
%--------------------------------------------------------------------------

% HP filter
populationraw.data = hpfilter(populationraw.data,10000);

%--------------------------------------------------------------------------
% Convert Monthly to Quarterly 
%--------------------------------------------------------------------------

Oilproductionraw = timeseries(mean(reshape(Oilproductionraw.data(1:end-2),3,[]))',datenum(datetime(1920,01,15):calmonths(3):datetime(2023,09,15))');

%--------------------------------------------------------------------------
% Synchronize Length
%--------------------------------------------------------------------------

[gdps,timespans]           = synchronize(gdpraw,timespan,'intersection');
[gdpdeflators,~]           = synchronize(gdpdeflatorraw,timespan,'intersection');
[consumptionnondurables,~] = synchronize(consumptionnondurableraw,timespan,'intersection');
[consumptionservicess,~]   = synchronize(consumptionservicesraw,timespan,'intersection');
[investments,~]            = synchronize(investmentraw,timespan,'intersection');
[investmentdeflators,~]    = synchronize(investmentdeflatorraw,timespan,'intersection');
[consumptiondurables,~]    = synchronize(consumptiondurableraw,timespan,'intersection');
[hourss,~]                 = synchronize(hoursraw,timespan,'intersection');
[fedfundsrates,~]          = synchronize(fedfundsrateraw,timespan,'intersection');
[populations,~]            = synchronize(populationraw,timespan,'intersection');
[HNcredits,~]              = synchronize(HNcreditraw,timespan,'intersection');
[HNrates,~]                = synchronize(HNrateraw,timespan,'intersection');
[NFcredits,~]              = synchronize(NFcreditraw,timespan,'intersection');
[NFrates,~]                = synchronize(NFrateraw,timespan,'intersection');
[oilprice,~]               = synchronize (oilpriceraw,timespan,'intersection');
[oil,~]                    = synchronize(oilproductionraw,timespan,'intersection');

%--------------------------------------------------------------------------
% Construct and Deflate Series
%--------------------------------------------------------------------------

gdp              = gdps.data ./ gdpdeflators.data;
consumption      = (consumptionnondurables.data + consumptionservicess.data) ./ gdpdeflators.data;
investment       = (investments.data + consumptiondurables.data) ./ gdpdeflators.data;
pinvest          = investmentdeflators.data ./ gdpdeflators.data;
HNcredit         = HNcredits.data ./ gdpdeflators.data;
NFcredit         = NFcredits.data ./ gdpdeflators.data;
HNspread         = HNrates.data - fedfundsrates.data;
NFspread         = NFrates.data - fedfundsrates.data;
Oilprice         = Oilprice.data ./ gdpdeflators.data; 

%--------------------------------------------------------------------------
% Express in Per Capita Terms
%--------------------------------------------------------------------------

gdp              = gdp              ./ populations.data;
consumption      = consumption      ./ populations.data;
investment       = investment       ./ populations.data;
hours            = hourss.data      ./ populations.data;
HNcredit         = HNcredit         ./ populations.data;
NFcredit         = NFcredit         ./ populations.data;
oilproduction    = Oilproduction    ./ pouplations.data;

%--------------------------------------------------------------------------
% Express in Log Difference and Remove Mean
%--------------------------------------------------------------------------

gdpgrowth         = 1 + log(gdp(2:end) ./ gdp(1:end-1)) - mean(log(gdp(2:end) ./ gdp(1:end-1)));
consumptiongrowth = 1 + log(consumption(2:end) ./ consumption(1:end-1)) - mean(log(consumption(2:end) ./ consumption(1:end-1)));
investmentgrowth  = 1 + log(investment(2:end) ./ investment(1:end-1)) - mean(log(investment(2:end) ./ investment(1:end-1)));
inflation         = 1 + log(gdpdeflators.data(2:end) ./ gdpdeflators.data(1:end-1)) - mean(log(gdpdeflators.data(2:end) ./ gdpdeflators.data(1:end-1)));
pinvestgrowth     = 1 + log(pinvest(2:end) ./ pinvest(1:end-1)) - mean(log(pinvest(2:end) ./ pinvest(1:end-1))); pinvesttrend = mean(log(pinvest(2:end) ./ pinvest(1:end-1)));
HNcreditgrowth    = 1 + log(HNcredit(2:end) ./ HNcredit(1:end-1)) - mean(log(HNcredit(2:end) ./ HNcredit(1:end-1)));
NFcreditgrowth    = 1 + log(NFcredit(2:end) ./ NFcredit(1:end-1)) - mean(log(NFcredit(2:end) ./ NFcredit(1:end-1)));
Oilproductiongrowth = 1 + log(Oilproduction(2:end) ./ Oilproduction(1:end-1)) - mean(log(Oilproduction(2:end) ./ Oilproduction(1:end-1)));

%--------------------------------------------------------------------------
% Express in Log and Remove Mean
%--------------------------------------------------------------------------

hours = 1 + log(hours(2:end)) - mean(log(hours(2:end)));
oilprice = 1 + log(oilprice(2:end)) - mean(log(oilprice(2:end)));

%--------------------------------------------------------------------------
% Express in Quarterly Percentage and Remove Mean
%--------------------------------------------------------------------------

fedfundsrate = 1 + (1 + fedfundsrates.data(2:end) / 100).^(1/4) - mean((1 + fedfundsrates.data(2:end) / 100).^(1/4));  
HNspread     = 1 + (1 + HNspread(2:end) / 100).^(1/4) - mean((1 + HNspread(2:end) / 100).^(1/4));
NFspread     = 1 + (1 + NFspread(2:end) / 100).^(1/4) - mean((1 + NFspread(2:end) / 100).^(1/4));


%--------------------------------------------------------------------------
% Rename Series
%--------------------------------------------------------------------------

gdp_obs              = gdpgrowth;
consumption_obs      = consumptiongrowth;
investment_obs       = investmentgrowth;
hours_obs            = hours;
inflation_obs        = inflation;
fedfundsrate_obs     = fedfundsrate;
pinvest_obs          = pinvestgrowth;
householdcredit_obs  = HNcreditgrowth; 
householdspread_obs  = HNspread;
firmcredit_obs       = NFcreditgrowth;
firmspread_obs       = NFspread;       
oil_obs              = Oilproductiongrowth;      
oilprice_obs         = oilprice;

%--------------------------------------------------------------------------
% Save Series
%--------------------------------------------------------------------------

% Levels
save('CSdatalevel','gdp','consumption','investment','hourss','gdpdeflators',...
     'fedfundsrates','pinvesttrend','HNcredit','NFcredit','HNspread','NFspread',...
      'oil','oilprice')

% Estimation
save('CSdataestim','gdp_obs','consumption_obs','investment_obs',... 
     'hours_obs','inflation_obs','fedfundsrate_obs','pinvest_obs',...
     'householdcredit_obs','householdspread_obs','firmcredit_obs',...
     'firmspread_obs','oil_obs','oilprice_obs')
