%% Oil supply, demand and risk shocks of Ready (2018), monthly, 1986 to date
%
% Ready, R. C. (2018). Oil prices and the stock market. Review of Finance
% 22(1), 155-176, splits each month's change in the oil price into three
% orthogonal parts using asset prices:
%   risk shock    v   the part explained by the unexpected change in the VIX
%   demand shock  d   the part explained by the return on oil producers'
%                     shares once the VIX innovation is taken out
%   supply shock  s   the rest
% With X = (oil price change, producer return, VIX innovation)' and
% Z = (s, d, v)', the paper sets X = A Z, A = [1 1 1; 0 a22 a23; 0 0 a33],
% and chooses A so that the shocks are uncorrelated (Section 3.1, equations
% 13-14). In every month s + d + v = the oil price change.
%
% Data (all free, downloaded when the script runs)
%   Oil price        WTI spot, month-end (FRED DCOILWTICO). The paper used the
%                    second-nearest NYMEX contract, which the EIA stopped
%                    publishing in April 2024; over 1986-2024 the two give
%                    supply shocks with a correlation of 0.97 and identical
%                    demand and risk shocks.
%   Producer return  Fama-French 49-industry "Oil" portfolio, value-weighted
%                    (Kenneth French's data library). The paper used the
%                    Datastream World Integrated Oil & Gas Producers index.
%   VIX              CBOE VIX, month-end (FRED VIXCLS), with the older VXO
%                    (FRED VXOCLS) for 1986-1989, before the VIX itself starts.
%   Recessions       NBER recession months (FRED USREC), for the chart only.
% Returns and price changes are logs; the VIX innovation is the residual of
% an ARMA(1,1) fitted to the log VIX over the whole sample, as in the paper.
%
% Output: oil_shocks_monthly.csv in the current folder, the same columns as
% update_oilshocks.py: month, oil_price, dp (log change in the oil price),
% xi_vix (VIX innovation), s, d, v (the shocks), dp_cum, s_cum, d_cum,
% v_cum (cumulated from January 1986), s_realtime, d_realtime, v_realtime
% (the shocks when the ARMA and the covariance matrix are estimated only on
% data up to each month, from a first window of 120 months), recession.
% The full-sample shocks are re-estimated every month, so earlier values
% move a little as the sample grows; the realtime columns never change.
%
% Cite: Ready, R. C. (2018). Oil prices and the stock market. Review of
% Finance 22(1), 155-176. https://doi.org/10.1093/rof/rfw071
%
% Requires MATLAB R2019b or later. No toolboxes: the ARMA(1,1) is estimated
% by conditional least squares with fminsearch, as in the Python script.

clear; clc;

START = datetime(1986,1,1);      % first month with a VIX-type index
VIX_FROM = datetime(1990,1,1);   % VIX from here, VXO before
FIRST_WINDOW = 120;              % months in the first realtime estimate

[tV, vix] = fred_month_end('VIXCLS');
[tO, vxo] = fred_month_end('VXOCLS');
[tW, wti] = fred_month_end('DCOILWTICO');
[tR, rec] = fred_monthly('USREC');
[tP, prod] = french_monthly('49_Industry_Portfolios_CSV.zip', 'Oil');

% Last month with every input
lastMonth = min([max(tV(~isnan(vix))), max(tW(~isnan(wti))), max(tP(~isnan(prod)))]);
month = (START:calmonths(1):lastMonth)';
n = numel(month);

VOL = align(tV, vix, month);
VXO = align(tO, vxo, month);
VOL(month < VIX_FROM) = VXO(month < VIX_FROM);
WTI = align(tW, wti, month);
PROD = align(tP, prod, month);
REC = align(tR, rec, month);

dp = [NaN; diff(log(WTI))];
rprod = log(1 + PROD);
lvix = log(VOL);
gaps = month(any(isnan([dp rprod lvix]), 2));
gaps = gaps(2:end);                                  % the first month has no price change
if ~isempty(gaps)
    error('Missing inputs in %s', strjoin(string(datestr(gaps, 'yyyy-mm')), ', '));
end

xi = arma11_residuals(lvix);
Z = nan(n, 3);
Z(2:end, :) = decompose([dp(2:end) rprod(2:end) xi(2:end)]);
s = Z(:,1); d = Z(:,2); v = Z(:,3);
cum = @(x) cumsum(fillmissing(x, 'constant', 0));

% Realtime: everything re-estimated on data up to each month
Zrt = nan(n, 3);
for k = FIRST_WINDOW:n
    xik = arma11_residuals(lvix(1:k));
    Zk = decompose([dp(2:k) rprod(2:k) xik(2:k)]);
    Zrt(k, :) = Zk(end, :);
end
xi(1) = NaN;

sig = @(x) round(x, 6, 'significant');
T = table(string(datestr(month, 'yyyy-mm')), sig(WTI), sig(dp), sig(xi), sig(s), sig(d), sig(v), ...
    sig(cum(dp)), sig(cum(s)), sig(cum(d)), sig(cum(v)), sig(Zrt(:,1)), sig(Zrt(:,2)), sig(Zrt(:,3)), REC, ...
    'VariableNames', {'month','oil_price','dp','xi_vix','s','d','v','dp_cum','s_cum','d_cum','v_cum', ...
    's_realtime','d_realtime','v_realtime','recession'});
writetable(T, 'oil_shocks_monthly.csv');
% Blank cells for missing values, as the Python script writes them
txt = fileread('oil_shocks_monthly.csv');
txt = regexprep(txt, '(?<=,)NaN(?=,|\r?\n|$)', '');
fid = fopen('oil_shocks_monthly.csv', 'w'); fwrite(fid, txt); fclose(fid);

share = var([s d v], 'omitnan') / var(dp, 'omitnan');
fprintf('Wrote oil_shocks_monthly.csv (%d months, %s to %s)\n', n, T.month(1), T.month(end));
fprintf('Share of the variance of oil price changes: supply %.0f%%, demand %.0f%%, risk %.0f%%\n', 100*share);
fprintf('%s: oil price change %+.1f%% = supply %+.1f%% + demand %+.1f%% + risk %+.1f%%\n', ...
    T.month(end), 100*dp(end), 100*s(end), 100*d(end), 100*v(end));

% Chart: cumulative shocks with NBER recession shading
figure('Color','w'); hold on;
C = 100 * [cum(s) cum(d) cum(v)];
yl = [floor(min(C(:))/50)*50, ceil(max(C(:))/50)*50];
inRec = REC == 1;
starts = find(diff([0; inRec]) == 1);  stops = find(diff([inRec; 0]) == -1);
for k = 1:numel(starts)
    x0 = month(starts(k));  x1 = month(stops(k)) + calmonths(1);
    patch([x0 x1 x1 x0], [yl(1) yl(1) yl(2) yl(2)], [0.85 0.85 0.85], 'EdgeColor','none', 'HandleVisibility','off');
end
plot(month, C(:,1), 'k', 'LineWidth', 1.4);
plot(month, C(:,2), 'Color', [0.14 0.25 0.56], 'LineWidth', 1.4);
plot(month, C(:,3), 'Color', [0.6 0.6 0.6], 'LineWidth', 1.2);
ylim(yl); grid on; box on;
ylabel('Cumulative shock since January 1986 (log points)');
legend({'Supply', 'Demand', 'Risk'}, 'Location', 'northwest');
title(sprintf('Oil shocks (Ready 2018), %s to %s', T.month(1), T.month(end)));

%% ---------------------------------------------------------------- method
function e = arma11_residuals(y)
% Residuals of y_t = c + phi y_{t-1} + theta e_{t-1} + e_t, estimated by
% conditional least squares (e_1 = 0). Same algorithm as the Python script.
    y = y(:);
    opts = optimset('TolX', 1e-8, 'TolFun', 1e-10, 'MaxFunEvals', 20000, 'MaxIter', 20000, 'Display', 'off');
    p = fminsearch(@(p) ssr(p, y), [0.1*mean(y) 0.9 0], opts);
    opts = optimset(opts, 'TolX', 1e-9, 'TolFun', 1e-12);
    p = fminsearch(@(p) ssr(p, y), p, opts);
    e = resid(p, y);
end

function e = resid(p, y)
    e = zeros(size(y));
    for t = 2:numel(y)
        e(t) = y(t) - p(1) - p(2)*y(t-1) - p(3)*e(t-1);
    end
end

function f = ssr(p, y)
    if abs(p(3)) >= 1, f = 1e12; return; end
    e = resid(p, y);
    f = sum(e(2:end).^2);
end

function Z = decompose(X)
% Columns of X: oil price change, producer return, VIX innovation. Returns
% Z = inv(A) X with A = [1 1 1; 0 a22 a23; 0 0 a33] chosen so that
% inv(A) cov(X) inv(A)' is diagonal (Ready 2018, equations 13-14). The
% producer return is demeaned first, so the demand and risk shocks average
% zero and the supply shock, the remainder, carries the average drift in
% the oil price; the covariances, and so A, are unaffected.
    S = cov(X);
    a33 = S(3,3) / S(1,3);
    a23 = S(2,3) / S(1,3);
    c = S(1,2) - S(2,3)*S(1,3)/S(3,3);
    e = S(2,2) - S(2,3)^2/S(3,3);
    a22 = e / c;
    A = [1 1 1; 0 a22 a23; 0 0 a33];
    X(:,2) = X(:,2) - mean(X(:,2));
    Z = (A \ X')';
end

%% ---------------------------------------------------------------- data
function [t, x] = fred_month_end(id)
% Last available daily value in each month of a FRED series.
    D = fred_csv(id);
    t = dateshift(D{:,1}, 'start', 'month');
    x = D{:,2};
    keep = ~isnan(x);
    t = t(keep); x = x(keep);
    [t, ~, g] = unique(t);
    last = accumarray(g, (1:numel(x))', [], @max);
    x = x(last);
end

function [t, x] = fred_monthly(id)
% A FRED series that is already monthly (dated the first of the month).
    D = fred_csv(id);
    t = dateshift(D{:,1}, 'start', 'month');
    x = D{:,2};
end

function D = fred_csv(id)
% FRED's server stalls on MATLAB's default request (its user agent, gzip
% encoding and Connection: close), so ask for an uncompressed reply on a
% kept-alive connection.
    url = ['https://fred.stlouisfed.org/graph/fredgraph.csv?id=' id];
    file = [tempname '.csv'];
    opts = weboptions('Timeout', 120, 'UserAgent', 'oilshocks-update/1.0', ...
        'HeaderFields', {'Accept-Encoding', 'identity'; 'Connection', 'keep-alive'});
    websave(file, url, opts);
    io = detectImportOptions(file);
    io = setvartype(io, 1, 'datetime');
    io = setvaropts(io, 1, 'InputFormat', 'yyyy-MM-dd');
    io = setvartype(io, 2, 'double');
    io = setvaropts(io, 2, 'TreatAsMissing', '.');
    D = readtable(file, io);
    delete(file);
end

function [t, x] = french_monthly(zipName, column)
% One column of the first monthly block (value-weighted returns, in percent)
% of a Kenneth French CSV file, as a decimal return.
    url = ['https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/' zipName];
    zf = [tempname '.zip'];
    websave(zf, url, weboptions('Timeout', 120, 'UserAgent', 'Mozilla/5.0'));
    folder = tempname; mkdir(folder);
    files = unzip(zf, folder);
    txt = fileread(files{1});
    delete(zf);
    lines = regexp(txt, '\r?\n', 'split');
    t = datetime.empty(0,1); x = [];
    col = 0;
    for i = 1:numel(lines)
        ln = lines{i};
        if col == 0
            if startsWith(ln, ',')
                head = strtrim(strsplit(ln, ',', 'CollapseDelimiters', false));
                c = find(strcmp(head, column), 1);
                if ~isempty(c), col = c; end
            end
            continue
        end
        tok = regexp(ln, '^\s*(\d{4})(\d{2}),', 'tokens', 'once');
        if isempty(tok), break; end
        cells = strsplit(ln, ',', 'CollapseDelimiters', false);
        t(end+1,1) = datetime(str2double(tok{1}), str2double(tok{2}), 1); %#ok<AGROW>
        x(end+1,1) = str2double(cells{col}); %#ok<AGROW>
    end
    x(x <= -99) = NaN;      % -99.99 marks missing
    x = x / 100;
end

function y = align(t, x, grid)
% Put a monthly series onto the monthly grid by exact month match.
    y = nan(size(grid));
    [hit, loc] = ismember(t, grid);
    y(loc(hit)) = x(hit);
end