backtest strategies using deep learning -凯发k8网页登录

pricedata = load('energyprices.mat','energyprices');
pricedata = pricedata.energyprices;
tail(pricedata)

serieslifespanplot(pricedata)

pricedata = removevars(pricedata,["nydiesel","gulfdiesel","laregular"]);

pricedata = rmmissing(pricedata);

returndata = tick2ret(pricedata)

% model is trained using a 30-day rolling window to predict 1 day in the
% future.
historysize = 30;
futuresize = 1;
% model predicts returns for oil, natural gas, propane, and kerosene.
outputvarname = ["brent" "naturalgas", "propane" "kerosene"];
numoutputs = numel(outputvarname);
% start_idx and end_idx are the index positions in the returndata
% timetable corresponding to the first and last date for making a prediction.
start_idx = historysize   1;
end_idx   = height(returndata) - futuresize   1;
numsamples = end_idx - start_idx   1;
% the date_vector variable stores the dates for making predictions.
date_vector = returndata.time(start_idx-1:end_idx-1);

network_features  = cell(numsamples,1);
network_responses = zeros(numsamples,numoutputs);
for j = 1:numsamples
    network_features{j} = (returndata(j:j historysize-1,:).variables)';
    network_responses(j,:) = ...
        (returndata(j historysize:j historysize futuresize-1,outputvarname).variables)';
end

% specify rows to use in the backtest (31-dec-2015 to 2-nov-2021).
backtest_start_idx = find(date_vector < datetime(2016,1,1),1,'last');
backtest_indices = backtest_start_idx:size(network_responses,1);
% specify data reserved for the backtest.
xbacktest = network_features(backtest_indices);
tbacktest = network_responses(backtest_indices,:);
% remove the backtest data.
network_features = network_features(1:backtest_indices(1)-1);
network_responses = network_responses(1:backtest_indices(1)-1,:);
% partition the remaining data into training and validation sets.
rng('default');
cv_partition = cvpartition(size(network_features,1),'holdout',0.2);
% training set
xtraining = network_features(~cv_partition.test,:);
ttraining = network_responses(~cv_partition.test,:);
% validation set
xvalidation = network_features(cv_partition.test,:);
tvalidation = network_responses(cv_partition.test,:);

numfeatures = width(returndata);
numhiddenunits_lstm = 10;
layers_lstm = [ ...
    sequenceinputlayer(numfeatures) 
    lstmlayer(numhiddenunits_lstm)
    layernormalizationlayer
    lstmlayer(numhiddenunits_lstm)
    layernormalizationlayer
    lstmlayer(numhiddenunits_lstm,'outputmode','last')
    layernormalizationlayer
    fullyconnectedlayer(numoutputs)
    regressionlayer];

max_epochs = 500;
mini_batch_size = 128;
learning_rate = 1e-4;
options_lstm = trainingoptions('adam', ...
    'plots','training-progress', ...
    'verbose',0, ...
    'maxepochs',max_epochs, ...
    'minibatchsize',mini_batch_size, ...
    'shuffle','every-epoch', ...
    'validationdata',{xvalidation,tvalidation}, ...
    'validationfrequency',50, ...
    'validationpatience',10, ...
    'initiallearnrate',learning_rate, ...
    'gradientthreshold',1);

dotrain = false;
if dotrain
    % train the lstm network.
    net_lstm = trainnetwork(xtraining,ttraining,layers_lstm,options_lstm);
else
    % load the pretrained network.
    load lstmbacktestnetwork
end

% compare the actual returns to model predicted returns.
actual    = tvalidation;
predicted = predict(net_lstm,xvalidation,'minibatchsize',mini_batch_size);

% overlay histogram of actual vs. predicted returns for the validation set.
output_idx = 1;
figure;
[~,edges] = histcounts(actual(:,output_idx),100);
histogram(actual(:,output_idx),edges);
hold on
histogram(predicted(:,output_idx),edges)
hold off
xlabel('percentage change in closing price')
legend('actual','predicted')
title(sprintf('%s: distribution of returns, actual vs. predicted', outputvarname(output_idx)))

% display the predicted vs. actual daily returns for the validation set.
figure
plot(actual(:,output_idx))
hold on
plot(predicted(:,output_idx))
yline(0)
legend({'actual','predicted'})
title(sprintf('%s: daily returns, actual vs. predicted', outputvarname(output_idx)))

% examine the residuals.
residuals = actual(:,output_idx) - predicted(:,output_idx);
figure;
normplot(residuals);

backtestpred_lstm = predict(net_lstm,xbacktest,'minibatchsize',mini_batch_size);

backtestsignaltt = timetable(date_vector(backtest_indices),backtestpred_lstm);

backtestpricett = pricedata(date_vector(backtest_indices),outputvarname);

risk_free_rate = 0.01;

% specify 10 basis points as the trading cost.
tradingcosts = 0.001;
% invest in long positions proportionally to their predicted return.
longstrategy = backteststrategy('longonly',@longonlyrebalancefcn, ...
    'transactioncosts',tradingcosts, ...
    'lookbackwindow',1);
% invest in both long and short positions proportionally to their predicted returns.
longshortstrategy = backteststrategy('longshort',@longshortrebalancefcn, ...
    'transactioncosts',tradingcosts, ...
    'lookbackwindow',1);
% invest 100% of capital into single asset with highest predicted returns.
bestbetstrategy = backteststrategy('bestbet',@bestbetrebalancefcn, ...
    'transactioncosts',tradingcosts, ...
    'lookbackwindow',1);
% for comparison, invest in an equal-weighted (buy low and sell high) strategy.
equalweightfcn = @(current_weights,prices,signal) ones(size(current_weights)) / numel(current_weights);
equalweightstrategy = backteststrategy('equalweight',equalweightfcn, ...
    'transactioncosts',tradingcosts, ...
    'lookbackwindow',0);

strategies = [longstrategy longshortstrategy bestbetstrategy equalweightstrategy];
bt = backtestengine(strategies,'riskfreerate',risk_free_rate);

bt = runbacktest(bt,backtestpricett,backtestsignaltt)

summary(bt)

figure;
equitycurve(bt)

function new_weights = longonlyrebalancefcn(current_weights,pricestt,signaltt) %#ok 
% long only strategy, in proportion to the signal.
signal = signaltt.backtestpred_lstm(end,:);
if any(0 < signal)
    signal(signal < 0) = 0;
    new_weights = signal / sum(signal);
else
    new_weights = zeros(size(current_weights));
end
end
function new_weights = longshortrebalancefcn(current_weights,pricestt,signaltt) %#ok 
% long/short strategy, in proportion to the signal
signal = signaltt.backtestpred_lstm(end,:);
abssum = sum(abs(signal));
if 0 < abssum
    new_weights = signal / abssum;
else
    new_weights = zeros(size(current_weights));
end
end
function new_weights = bestbetrebalancefcn(current_weights,pricestt,signaltt) %#ok 
% best bet strategy, invest in the asset with the most upside.
signal = signaltt.backtestpred_lstm(end,:);
new_weights = zeros(size(current_weights));
new_weights(signal == max(signal)) = 1;
end
function serieslifespanplot(pricedata)
% plot the lifespan of each time series.
% specify all time series end on same day.
d2 = numel(pricedata.time);
% plot the lifespan patch for each series.
numseries = size(pricedata,2);
for i = 1:numseries
    % find start date index.
    d1 = find(~isnan(pricedata{:,i}),1,'first');
    % plot patch.
    x = [d1 d1 d2 d2];
    y = i   [-0.4 0.4 0.4 -0.4];
    patch(x,y,[0 0.4470 0.7410])
    hold on
end
hold off
% set the plot properties.
xlim([-100 d2]);
ylim([0.2 numseries   0.8]);
yticks(1:numseries);
yticklabels(pricedata.properties.variablenames');
flipud(gca);
years = 1990:5:2021;
xtick_idx = zeros(size(years));
for yidx = 1:numel(years)
    xtick_idx(yidx) = find(years(yidx) == year(pricedata.time),1,'first');
end
xticks(xtick_idx);
xticklabels(string(years));
title('time series life span');
end