Skip to content

Commit

Permalink
Split off generating uncertainty into dedicated function.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lbrandt
lbrandt committed Feb 16, 2018
1 parent 2c0f7ef commit 1c8457d
Show file tree
Hide file tree
Showing 2 changed files with 214 additions and 132 deletions.
152 changes: 152 additions & 0 deletions MATLAB/aggregateUncertainty.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
function results = generateUncertainty(sy, ty, sf, tf, ybeta, py, fbeta, pf, h)
% -------------------------------------------------------------------------
% Generates aggregate uncertainty embodied in a set of prediction errors
% from a factor model together with their latent stochastic volatility.
% Adapted from Jurado, Ludvigson, Ng (2015).
%
% Input
% sy Stochastic volatility in macro variables y [T x N]
% ty AR model parameters of sy process [3 x N]
% sf Stochastic volatility in factors y [T x R]
% tf AR model parameters sf process [3 x R]
% ybeta Parameters of prediction models for y [I x N]
% py Number of autoregressive lags in ybeta
% fbeta Parameters of prediction models for f [J x N]
% pf Number of autoregressive lags in fbeta
% h Maximum number of forecasting steps
%
% Output
% results Structure containing:
% results.uavg Aggregate uncertainty by simple average [T x h]
% results.ufac Aggregate uncertainty by PCA weighting [T x h]
% results.uind Estimates of individual expected volatilities [T x N x h]
% results.phic Cell array of parameter matrices {N}
%
% Dependencies {source}
% companion
% expectvar
%
% -------------------------------------------------------------------------

%%%%
% Initialisation
T = length(sy);

bf = sparse(fbeta);
by = sparse(ybeta);

[numparf, R] = size(bf);
[numpary, N] = size(by);

% Reparameterise means
tf(1, :) = tf(1, :).* (1 - tf(2, :));
ty(1, :) = ty(1, :).* (1 - ty(2, :));


%%%%
% Compute expected volatility in factors

% Expected h-step-ahead variance Et[sigma2_F(t+h)]
evarf = zeros(T, R, h);
for j = 1:h
for i = 1:R
alpha = tf(1,i);
beta = tf(2,i);
tau = tf(3,i);
x = sf(:,i);

evarf(:, i, j) = expectvar(x, alpha, beta, tau, j);
end
end


% Build VAR representation of the system of factors
fvar = bf(1, :)'; % Collect intercepts
for i = 2:numparf % Append VAR parameter matrices, they are diagonal because factors are not cross-correlated
fvar = [fvar, diag(bf(i, :))]; %#ok<AGROW> % Suppress preallocation error
end

% Build parameter matrix Phi of companion form
phif = companion(R, pf, fvar);


%%%%
% Compute uncertainty in macro variables

evary = zeros(T, N, h);
ut = zeros(T, N, h);
results.phic = cell(N, 1);
for i = 1:N

tic;

% Build parameter matrix for variable's FAVAR representation
lambda = sparse(1, 1:numpary-(py+1), by(py+2:end, i), py, R*pf); % Lambda contains parameters from reg(y ~ z) on first row, then filled up with zeros to match dims
phiy = companion(1, py, by(1:py+1, i)');

phi = [phif, zeros(R*py, py); lambda, phiy];


% Compute individual uncertainty via recursion
alpha = ty(1, i);
beta = ty(2, i);
tau = ty(3, i);
x = sy(:, i);


% Compute h-step-ahead expected variance in variable i Et[sigma2_Y(t+h)]
for j = 1:h
evary(:, i, j) = expectvar(x, alpha, beta, tau, j);
end

for t = 1:T
for j = 1:h

% Diagonal matrix of h-step-ahead variance forecasts for each variable/factor separately
% Zeros on diag where depvar vector [Zt, Yjt]' contains lagged terms
evh = sparse(1:R*pf + py, 1:R*pf + py, [evarf(t, :, j)'; zeros(R*pf - R, 1); evary(t, i, j)'; zeros(py - 1, 1)]);

if j == 1
ui = evh; % for h=1, Omega = Et[sigma2_y(t+h)]
else
ui = phi* ui* phi' + evh; % for h>1, Omega = phi* Et[sigma2_Y(t+h-1)]* phi + Et[sigma2_Y(t+h)]
end

results.phic{i} = phi; % save phi matrices

ut(t, i, j) = ui(R*pf + 1, R*pf + 1); % Select element corresponding to yt, this is the squared uncertainty in this variable at time t looking h-steps ahead
end
end

fprintf('Series %d, elapsed time = %0.4f \n', i, toc);

end
% Individual expected volatilities
results.uind = sqrt(ut);



%%%%
% Aggregate individual uncertainty
% Simple average
results.uavg = squeeze(mean(results.uind,2));

% Principal component analysis on logged variances
logu = log(ut(:, :, :));

results.ufac = zeros(T, h);
upca = zeros(T, h);
for i = 1:h
upca(:, i) = factors(logu(:, :, i), 1, 2);

% Flip ufac if necessary
rho = corrcoef(upca(:, i), results.uavg(:, i));
if rho(2, 1) < 0
upca = -upca;
end

% Scale to Uavg
results.ufac(:, i) = exp( standardise(upca(:, i))* std(log(results.uavg(:, i))) + mean(log(results.uavg(:, i))) );
end

end
194 changes: 62 additions & 132 deletions MATLAB/ez_uncertainty.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,157 +8,87 @@
% Load data
load ez_factors_forc;

%h5disp('gs_svfresults.h5')
sf = h5read('ez_svfresults.h5', '/h')'; % Estimates of latent process s(t)
tf = h5read('ez_svfresults.h5', '/t')'; % AR(1) parameter estimators
% New
%h5disp('ez_svresults.h5')
sf = h5read('ez_svresults.h5', '/sf')'; % Estimates of latent process s(t)
tf = h5read('ez_svresults.h5', '/tf')'; % AR(1) parameter estimators

%h5disp('gs_svyresults.h5')
sy = h5read('ez_svyresults.h5', '/h')';
ty = h5read('ez_svyresults.h5', '/t')';
% From LASSO model selection
sy = h5read('ez_svresults.h5', '/sy')';
ty = h5read('ez_svresults.h5', '/ty')';


%%%%
% Compute uncertainty in macro variables


% Initialisation
h = 12;

T = size(sy, 1);

bf = sparse(fbetas);
by = sparse(ybetas);

[numparf, R] = size(bf);
[numpary, N] = size(by);


tf(1, :) = tf(1, :).* (1 - tf(2, :)); % Reparameterise mean
ty(1, :) = ty(1, :).* (1 - ty(2, :)); % Reparameterise mean
% From hard thresholding model selection
syht = h5read('ez_svresults.h5', '/syht')';
tyht = h5read('ez_svresults.h5', '/tyht')';



%%%%
% Compute expected volatility in factors

% Compute uncertainty in macro variables
hmax = 12;

% Expected h-step-ahead variance Et[sigma2_F(t+h)]
evarf = zeros(T, R, h);
uLasso = aggregateUncertainty(sy, ty, sf, tf, ybetas, py, fbetas, pf, hmax);
uHT = aggregateUncertainty(syht, tyht, sf, tf, htbetas, py, fbetas, pf, hmax);

for j = 1:h
for i = 1:R
alpha = tf(1,i);
beta = tf(2,i);
tau = tf(3,i);
x = sf(:,i);

evarf(:, i, j) = expectvar(x, alpha, beta, tau, j);
end
end

% Compare methods

% Build VAR representation of the system of factors
% Plots
figure
for i = 1:hmax
subplot(4, 3, i)
plot(dates, uLasso.uavg(:, i))
hold on
plot(dates, uHT.uavg(:, i))
title(['h = ',num2str(i)])
end

fvar = bf(1, :)'; % Collect intercepts
for i = 2:numparf % Append VAR parameter matrices, they are diagonal because factors are not cross-correlated

fvar = [fvar, diag(bf(i, :))]; %#ok<AGROW> % Suppress preallocation error
figure
for i = 1:hmax
subplot(4, 3, i)
plot(dates, uLasso.ufac(:, i))
hold on
plot(dates, uHT.ufac(:, i))
title(['h = ',num2str(i)])
end

% Build parameter matrix Phi of companion form
phif = companion(R, pf, fvar);
% Means
ubar1 = mean(uLasso.uavg);
ubar2 = mean(uHT.uavg);

ubardiff = ubar1 - ubar2;

%%%%
% Compute uncertainty in macro variables


% Initialisation
evary = zeros(T, N, h);
ut = zeros(T, N, h);

phiy = sparse(py, py, 0);
phi = sparse(R*pf + py, R*pf + py, 0);

phisav = cell(N,1);

for i = 1:N

tic;

% Build parameter matrix for variable's FAVAR representation
lambda = sparse(1, 1:numpary-(py+1), by(py+2:end, i), py, R*pf); % Lambda contains parameters from reg(y ~ z) on first row, then filled up with zeros to match dims
phiy = companion(1, py, by(1:py+1, i)');

phi = [phif, zeros(R*py, py); lambda, phiy];


% Compute individual uncertainty via recursion
alpha = ty(1, i);
beta = ty(2, i);
tau = ty(3, i);
x = sy(:, i);


% Compute h-step-ahead expected variance in variable i Et[sigma2_Y(t+h)]
for j = 1:h
evary(:, i, j) = expectvar(x, alpha, beta, tau, j);
end

for t = 1:T
for j = 1:h

% Diagonal matrix of h-step-ahead variance forecasts for each variable/factor separately
% Zeros on diag where depvar vector [Zt, Yjt]' contains lagged terms
evh = sparse(1:R*pf + py, 1:R*pf + py, [evarf(t, :, j)'; zeros(R*pf - R, 1); evary(t, i, j)'; zeros(py - 1, 1)]);

if j == 1
ui = evh; % for h=1, Omega = Et[sigma2_y(t+h)]
else
ui = phi* ui* phi' + evh; % for h>1, Omega = phi* Et[sigma2_Y(t+h-1)]* phi + Et[sigma2_Y(t+h)]
end

phisav{i} = phi; % save phi matrices

ut(t, i, j) = ui(R*pf + 1, R*pf + 1); % Select element corresponding to yt, this is the squared uncertainty in this variable at time t looking h-steps ahead
end
end

fprintf('Series %d, Elapsed Time = %0.4f \n', i, toc);

%%%%
% Save results
selectMethod = uLasso;
uavg = selectMethod.uavg;
ufac = selectMethod.ufac;
uing = selectMethod.uind;
phic = selectMethod.phic;
%save ez_uncertainty -v7.3 dates uavg ufac uing phic


% Figure from JLN2015
selectH = [1, 3, 12];

figure
subplot(2,1,1);
for i = selectH
plot(dates, uavg(:, i), 'DisplayName', ['h = ',num2str(i)])
hold on
end
legend('show')
title('uavg')

% Individual expected volatilities
Uind = sqrt(ut);


%%%%
% Aggregate individual uncertainty
% Simple average
Uavg = squeeze(mean(Uind,2));


% Principal component analysis on logged variances
logut = log(ut(:, :, :));

Upca = zeros(T, h);
Ufac = zeros(T, h);

for i = 1:h

Upca(:, i) = factors(logut(:, :, i), 1, 2);

% Flip Ufac if necessary
rho = corrcoef(Upca(:, i), Uavg(:, i));
if rho(2, 1) < 0
Upca = -Upca;
end

% Scale to Uavg
Ufac(:, i) = exp( standardise(Upca(:, i))* std(log(Uavg(:, i))) + mean(log(Uavg(:, i))) );
subplot(2,1,2);
for i = selectH
plot(dates, ufac(:, i), 'DisplayName', ['h = ',num2str(i)])
hold on
end
legend('show')
title('ufac')



%%%%
% Save results
save ez_uncertainty dates Uavg Ufac ut phisav

0 comments on commit 1c8457d

Please sign in to comment.