Backtest a Profitable Trend-Following Strategy

In this article, we provide a detailed walkthrough of the Matlab code used to perform the trend-following backtest discussed in our paper, A Century of Profitable Industry Trends.

We recommend reading the paper first, which you can access here.

Using Kenneth French’s freely available database, we built an industry-based long-only trend-following portfolio. This portfolio is based on daily data from 48 industry portfolios, covering the period from 1926 to 2024. Our study compares the performance of this momentum-based strategy with a passive Buy & Hold approach over the past century.

Our findings demonstrate the superior performance of momentum-based portfolios. The strategy includes several parameters, none of which were optimized in-sample, meaning the performance statistics could be improved by fine-tuning these parameters.

At the end of this article, we provide another Matlab code that allows you to quickly explore how changing different parameters can impact the portfolio’s profitability.

To enhance the readability and make the code more intuitive for novice quant researchers, we intentionally avoided more complex coding procedures that may result in improved computational efficiency.

Step-by-Step Guide Through the Backtesting Process

Below, we give an overview of the main building blocks behind the backtesting procedure.

Step 1: Download and Process Industry Data

Step 2: Load Database and Compute Indicators and Signals

Step 3: Run the Backtest

Step 4: Study Results

Tools Needed

MATLAB: MATLAB is a high-level programming language and environment developed by MathWorks. It is widely used for numerical computing, data analysis, algorithm development, and visualization. You can try MATLAB for FREE for 20 hours by visiting the official MathWorks website. Alternatively, a MATLAB Home license costs $149 for personal use.

How to Read This Post

For each step, we will provide the full code first, then explain it step by step. If a step depends on a utility function (a function we wrote that is not built-in MATLAB), we will provide its code after the step code directly before the explanation.

Step 1: Download and Process Industry Data

Overview

This step involves downloading daily return data for 48 industry portfolios from Kenneth French‘s database, processing the data to prepare it for backtesting, and saving the necessary variables for later use. The process includes downloading, unzipping, reading CSV data, cleaning, and structuring the data.

Step 1 Code:

Click to see the MATLAB code for Step 1

MATLAB

% Download Data From French
url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily_CSV.zip';

zipFile = websave('temp.zip', url);
pause(1);

unzip(zipFile, 'temp_folder');

pause(1);

csvFile = dir(fullfile('temp_folder', '*.csv'));
csvFilePath = fullfile('temp_folder', csvFile.name);

movefile(csvFilePath, 'output.csv', 'f');

delete(zipFile); pause(1);
rmdir('temp_folder', 's');

clear p 

csvData          = readtable('output.csv');
p.industry_names = csvData.Properties.VariableNames(2:end);
rows_nan         = find(isnan(table2array(csvData(:,1))));

% MktCap ind_weighting
from     = 1; until = rows_nan(1)-1;
data     = table2array(csvData(from:until,1:end)); 
p.ret    = data(:,2:end)/100;
p.caldt  = datenum(num2str(data(:,1)),'yyyymmdd'); % change in case daily data
p.ret(find(p.ret<=-0.99))=NaN;

%%

clearvars -except p

p.mkt_ret   = market_french_reconciled(p.caldt); % download the time-series of market return and match the position with respect to days
p.tbill_ret = tbill_french_reconciled(p.caldt);  % download the time-series of tbill daily ret and match the position with respect to days

save('industry48.mat','p')

Needed functions for Step 1:

market_french_reconciled
tbill_french_reconciled

Click to see the MATLAB code for market_french_reconciled

MATLAB

function y = market_french_reconciled(caldt)

url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_CSV.zip';

zipFile = websave('temp.zip', url);
pause(1);

unzip(zipFile, 'temp_folder');

pause(1);

csvFile = dir(fullfile('temp_folder', '*.csv'));
csvFilePath = fullfile('temp_folder', csvFile.name);
pause(1);

movefile(csvFilePath, 'output.csv', 'f');

delete(zipFile); pause(2);
rmdir('temp_folder', 's');


csvData = readtable('output.csv', 'ReadVariableNames', false);
rows_nan = find(isnan(table2array(csvData(:,1))));

% MktCap Weighting
from     = rows_nan(1); until = rows_nan(2)-1;
data = table2array(csvData(from:until,1:end)); 

data(find(isnan(data(:,1))),:)=[];
ret_mkt        = data(:,2)/100+data(:,5)/100; % add back the risk free
caldt_mkt      = datenum(num2str(data(:,1)),'yyyymmdd'); % change in case daily data

T = length(caldt);

y = NaN(T,1);

X_dates = datetime(caldt, 'ConvertFrom', 'datenum');


[isMember, idx] = ismember(caldt_mkt, caldt);
% Create a mapping of caldt_mkt to ret_mkt using idx
mappedReturns = NaN(size(caldt));
mappedReturns(idx(isMember)) = ret_mkt(isMember);
% Assign the returns to y based on matching indices
y = mappedReturns;


end

Click to see the MATLAB code for tbill_french_reconciled

MATLAB


function y = tbill_french_reconciled(caldt)

url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_CSV.zip';

zipFile = websave('temp.zip', url);
pause(1);

unzip(zipFile, 'temp_folder');

pause(1);

csvFile = dir(fullfile('temp_folder', '*.csv'));
csvFilePath = fullfile('temp_folder', csvFile.name);
pause(1);

movefile(csvFilePath, 'output.csv', 'f');

delete(zipFile); pause(2);
rmdir('temp_folder', 's');


csvData = readtable('output.csv', 'ReadVariableNames', false);
rows_nan = find(isnan(table2array(csvData(:,1))));

% MktCap Weighting
from     = rows_nan(1); until = rows_nan(2)-1;
data = table2array(csvData(from:until,1:end)); 

data(find(isnan(data(:,1))),:)=[];
ret_mkt        = data(:,5)/100; % 
caldt_mkt      = datenum(num2str(data(:,1)),'yyyymmdd'); % change in case daily data

T = length(caldt);

y = NaN(T,1);

X_dates = datetime(caldt, 'ConvertFrom', 'datenum');




[isMember, idx] = ismember(caldt_mkt, caldt);
% Create a mapping of caldt_mkt to ret_mkt using idx
mappedReturns = NaN(size(caldt));
mappedReturns(idx(isMember)) = ret_mkt(isMember);
% Assign the returns to y based on matching indices
y = mappedReturns;

end

1.1 Download and Store Data

This code downloads a ZIP file from the specified URL, which contains CSV files for 48 industry portfolios. The websave function is used to download and save the file as ‘temp.zip’. The pause commands are used to ensure that file operations complete before proceeding. The ZIP file is then extracted into a folder named ‘temp_folder’.

Click to see the MATLAB code

MATLAB

url = 'https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily_CSV.zip';
zipFile = websave('temp.zip', url);
pause(1);
unzip(zipFile, 'temp_folder');
pause(1);

1.2 Locate and Organize Data Files

After extracting the files, the script searches for CSV files within the ‘temp_folder’ and moves the found CSV file to the current working directory, renaming it to ‘output.csv’. This ensures that subsequent operations are performed on the correct file.

Click to see the MATLAB code

MATLAB

csvFile = dir(fullfile('temp_folder', '*.csv'));
csvFilePath = fullfile('temp_folder', csvFile.name);
movefile(csvFilePath, 'output.csv', 'f');

1.3 Clean Up Temporary Files

This part of the script deletes the original ZIP file and removes the temporary directory created for extraction, cleaning up the workspace.

Click to see the MATLAB code

MATLAB

delete(zipFile); pause(1);
rmdir('temp_folder', 's');

1.4 Load Data and Initialize Variables

Clears variable p to ensure it does not contain old data. The script then reads the processed CSV file into a table, extracts industry names, and identifies rows containing NaN values which could indicate missing data.

Click to see the MATLAB code

MATLAB

clear p
csvData = readtable('output.csv');
p.industry_names = csvData.Properties.VariableNames(2:end);
rows_nan = find(isnan(table2array(csvData(:,1))));

1.5 Prepare Financial Data

Converts the table to an array for easier manipulation, slicing the data until the first NaN occurrence. It normalizes the returns by dividing by 100 to convert percentage to a decimal format. Dates are converted to MATLAB’s numeric date format for consistency in subsequent operations. Extreme negative returns, usually associated with missing data, are marked as NaN to avoid distorting statistical calculations.

Click to see the MATLAB code

MATLAB

from = 1; until = rows_nan(1)-1;
data = table2array(csvData(from:until,1:end)); 
p.ret = data(:,2:end)/100;
p.caldt = datenum(num2str(data(:,1)),'yyyymmdd');
p.ret(find(p.ret<=-0.99))=NaN;

1.6 Download and Align Additional Data

Calls custom functions to download and align market and T-bill returns with the dates from the industry data. These functions are crucial for ensuring that all data used in subsequent analysis is synchronized to the same time frame.

Click to see the MATLAB code

MATLAB

p.mkt_ret = market_french_reconciled(p.caldt);
p.tbill_ret = tbill_french_reconciled(p.caldt);

1.7 Save Processed Data

The final step saves the prepared dataset to a MATLAB file, ‘industry48.mat’. This file will be used in subsequent steps of the backtesting process.
Saving the database allows the quant researcher to import data quicker for future backtest.

Click to see the MATLAB code

MATLAB

save('industry48.mat','p')

Step 2: Load Database and Compute Indicators and Signals

Overview

In this step, after loading the pre-processed financial data (industry48.mat), various technical indicators and trading signals are computed. This includes setting up the parameters for these indicators, generating hypothetical price series from total returns, and calculating volatility, moving averages, Donchian channels, Keltner bands, and the corresponding long and short trading signals. These calculations are foundational for the strategy’s decision-making process in the subsequent backtesting phase.

Step 2 Code:

Click to see the MATLAB code

MATLAB

clear all

load('industry48.mat')

% INDICATORS PARAMETERS
up_day      = 20;
down_day    = 40;
kelt_mult   = 2;
adr_vol_adj = 1.4; % we noticed that ATR is usually 1.4xVol(close2close)
kelt_mult   = kelt_mult*adr_vol_adj;

% FROM DAILY TOTAL RETURNS, CREATE THE HYPOTHETICAL PRICE TIME-SERIES FOR
% EACH INDUSTRY
p.price       = cumprod(1+p.ret,'omitmissing');

% Rolling Volatility of Daily Returns
p.vol         = Rolling_Vol(p.ret,up_day,'omtinan');

% Technical Indicators
p.ema_down    = Rolling_Ema(p.price,down_day,'omitnan');
p.ema_up      = Rolling_Ema(p.price,up_day,'omitnan') ;

% Donchain Channels (Eq.5 and Eq.6)
p.donc_up     = Rolling_Max(p.price,up_day);
p.donc_down   = Rolling_Min(p.price,down_day);

% Keltner Bands (Eq. 3 and 4)
p.kelt_up     = p.ema_up + kelt_mult*Rolling_Mean(abs(price2change(p.price,1)),up_day,'omitnan');
p.kelt_down   = p.ema_down - kelt_mult*Rolling_Mean(abs(price2change(p.price,1)),down_day,'omitnan');

% Model Bands 
p.long_band   = min(p.donc_up,p.kelt_up); % Eq.7
p.short_band  = max(p.donc_down,p.kelt_down); % Eq.11

% Model Long Signal
p.long_signal = p.price>=lag_TS(p.long_band,1) & lag_TS(p.long_band,1)>lag_TS(p.short_band,1);

Needed functions for Step 2:

Rolling_Vol
Rolling_Ema
Rolling_Max
Rolling_Min
Rolling_Mean
price2change
lag_TS

Click to see the MATLAB code for Rolling_Vol

MATLAB

function y=Rolling_Vol(price,window,nanflag)
T = size(price,1);
N = size(price,2);

y = NaN(T,N);

if prod(double(nanflag=='omitnan'))==1
    
            for i=window:1:T
                y(i,:)=std(price(i-window+1:i,:),1,'omitmissing');
            end
else
            for i=window:1:T
                y(i,:)=std(price(i-window+1:i,:),1,'includenan');
            end

end        
end

Click to see the MATLAB code for Rolling_Ema

MATLAB

function ema = rolling_ema(price, window, nanflag)

% Size of the price matrix
T = size(price, 1);
N = size(price, 2);

% Initialize output matrix
ema = NaN(T, N);

% Smoothing factor
alpha = 2 / (window + 1);

% Loop through each element and compute EMA
for n = 1:N
    for t = 2:T
        if isnan(price(t, n))
            if strcmp(nanflag, 'omitnan')
                % If omitnan flag is used, skip NaN values
                continue;
            end
        else
            % If the previous EMA is NaN (start of the series), initialize it as the first data point
            if isnan(ema(t-1, n))
                ema(t-1, n) = price(t-1, n);
            end
            % Compute the EMA
            ema(t, n) = alpha * price(t, n) + (1 - alpha) * ema(t-1, n);
        end
    end
end

end

Click to see the MATLAB code for Rolling_Max

MATLAB

function y=Rolling_Max(price,window)
y = movmax(price,[window-1 0],1,"omitmissing");
if size(y,1) >= window
y(1:window-1,:) = NaN;
end

end

Click to see the MATLAB code for Rolling_Min

MATLAB

function y=Rolling_Min(price, window)

y = movmin(price,[window-1 0],1,"omitmissing");
if size(y,1) >= window
y(1:window-1,:) = NaN;
end

end

Click to see the MATLAB code for Rolling_Mean

MATLAB

function y=Rolling_Mean(price,window,nanflag)

T = size(price,1);
N = size(price,2);

y = NaN(T,N);

  
for i=window:1:T
    y(i,:)=mean(price(i-window+1:i,:),1,nanflag);
end

y(find(isnan(price)==1))=NaN; % if input variable is missing, consider NaN
        
end

Click to see the MATLAB code for price2change

MATLAB

function y=price2change(price,n)

    % This function converts prices into returns

    T=size(price,1);
    N=size(price,2);

    y=NaN(T,N);
    
    
    y(n+1:end,:)=price(n+1:T,:)-price(1:T-n,:);

    y(isnan(y))=NaN;    
    
end

Click to see the MATLAB code for lag_TS

MATLAB

function y=lag_TS(data,n)

    % This function lag the data matrix to n-times

    T=size(data,1);
    N=size(data,2);

    y=NaN(T,N);
    
    if n>0
     y(n+1:T,:)   = data(1:T-n,:);
    elseif n<0
        n = -n;
        y(1:T-n,:) = data(1+n:T,:);   
    else 
        y = data;
    end

    y(isnan(y))=NaN;    
    
end

2.1 Load Data

This section begins by clearing the MATLAB workspace with clear all to ensure no residual data from previous sessions affects the current process. It then loads the dataset ‘industry48.mat’ which was prepared in Step 1.

Click to see the MATLAB code

MATLAB

clear all
load('industry48.mat')

2.2 Set Indicator Parameters

Parameters for the technical indicators are defined here. up_day and down_day specify the periods for calculating moving averages and other metrics. The Keltner Channel multiplier (kelt_mult) is multiplied by an adjustment factor (adr_vol_adj) since the database we are using does not include high and low prices (more information can be found on page 7 of our paper)

Click to see the MATLAB code

MATLAB

up_day      = 20;
down_day    = 40;
kelt_mult   = 2;
adr_vol_adj = 1.4; % we noticed that ATR is usually 1.4xVol(close2close)
kelt_mult   = kelt_mult*adr_vol_adj;

2.3 Compute Price Series

Constructs a hypothetical price time-series from daily total returns for each industry by sequentially multiplying the compounded returns. This is essential for the subsequent calculation of technical indicators which require price data.

Click to see the MATLAB code

MATLAB

p.price = cumprod(1+p.ret,'omitmissing');

2.4 Calculate Volatility

Calculates rolling volatility over a specified number of days (up_day). This metric is crucial for risk management and sizing positions based on current market conditions.

Click to see the MATLAB code

MATLAB

p.vol = Rolling_Vol(p.ret, up_day, 'omitnan');

2.5 Generate Technical Indicators

This section computes various technical indicators:
- Exponential Moving Averages (EMA) for up and down days.
- Donchian Channels which are the maximum and minimum of the price over a specified period.
- Keltner Channels which use the EMA and the mean absolute price change to set channel boundaries.

Click to see the MATLAB code

MATLAB

p.ema_down = Rolling_Ema(p.price, down_day, 'omitnan');
p.ema_up = Rolling_Ema(p.price, up_day, 'omitnan');
p.donc_up = Rolling_Max(p.price, up_day);
p.donc_down = Rolling_Min(p.price, down_day);
p.kelt_up = p.ema_up + kelt_mult * Rolling_Mean(abs(price2change(p.price, 1)), up_day, 'omitnan');
p.kelt_down = p.ema_down - kelt_mult * Rolling_Mean(abs(price2change(p.price, 1)), down_day, 'omitnan');

2.6 Determine Trading Signals

Trading signals are derived by comparing the current price against the previously computed trading bands:
- Long Signal: Generated when the current closing price is above the long band as of yesterday’s close. To prevent contradicting signals, we also require that yesterday’s long band should be above yesterday’s short band.

Click to see the MATLAB code

MATLAB

p.long_band = min(p.donc_up, p.kelt_up); 
p.short_band = max(p.donc_down, p.kelt_down); 
p.long_signal = p.price >= lag_TS(p.long_band, 1) & lag_TS(p.long_band, 1) > lag_TS(p.short_band, 1);

Step 3: Run the Backtest

Overview

Step 3 involves executing the backtest using the previously calculated indicators and signals to simulate trading on a historical data set. It focuses on dynamic portfolio management based on the computed technical indicators and trading signals, with adjustments made for exposure, position sizing, and stop-loss levels. The process includes setting up initial parameters such as investment capital and risk measures, dynamically updating exposures and weights based on trading signals, and calculating returns to measure performance.

Step 3 Code

Click to see the MATLAB code

MATLAB

AUM_0           = 1;

invest_cash     = "YES"; % invest cash in short-term bonds, and pay for cash borrowed

target_vol      = 0.02; % the risk budget per industry is target_vol/N
max_leverage    = 2;    % max leverage at portfolio level.
max_not_trade   = 0.20; % Limit the max equity exposure per industry

% initiliaze portfolio variables
p.exposure         =  zeros(size(p.price)); % matrix that take 1 or 0 based on if we hold a long position
p.ind_weight       =  zeros(size(p.price)); % optimal weight at industry level
p.trail_stop_long  =  NaN(size(p.price)); % trailing stop loss

N_ind              = length(p.industry_names); % how many industries in database

% for each industry, compute the optimal weight at daily frequency, and
% update the level of trailing stop in case the stop has not been hit

T = length(p.caldt);

for j = 1:length(p.industry_names)

    for t = 1:T
        
        if ~isnan(p.ret(t,j)) && ~isnan(p.long_band(t,j)) % start the backtest when key-variable are not NaN

            % EXPOSURE CALCULATION
            
            % Open a new long position if yesterday exposure was 0 and
            % today we have a long_signal at the closure
            if p.exposure(t-1,j) <= 0 && p.long_signal(t,j) == 1 % New long signal
               
               p.exposure(t,j)        = 1;
               p.trail_stop_long(t,j) = p.short_band(t,j); % set the value for the trailing stop
                              
               % leverage required to reach the target vol 
               lev_vol            = target_vol/p.vol(t,j);
               
               % optimal weight at industry level
               p.ind_weight(t,j)  = p.exposure(t,j)*lev_vol;

               
            % Yesterday we were long. Today the price didn't cross the trailing stop. We confirm the long epxosure and we update the trailing stop   
            elseif p.exposure(t-1,j) == 1 && p.price(t,j) > max(p.trail_stop_long(t-1,j),p.short_band(t,j))  

               p.exposure(t,j)        = 1; % the long exposure is confirmed
               p.trail_stop_long(t,j) = max(p.trail_stop_long(t-1,j),p.short_band(t,j)); % update trailing stop as in Eq.12

               % leverage required to reach the target vol 
               lev_vol           = target_vol/p.vol(t,j);

               % optimal weight at industry level
               %p.ind_weight(t,j) = min(p.ind_weight(t-1,j),lev_vol); % in the paper i used this but it's a small typo, it deflated
               %the results in terms of IRR
               p.ind_weight(t,j) = p.exposure(t,j)*lev_vol;
               
            % yesterday we were long but today the closing price is below the trailing stop   
            elseif p.exposure(t-1,j) == 1 && p.price(t,j) <= max(p.trail_stop_long(t-1,j),p.short_band(t,j)) 
               
               % set exposure and weight to 0
               p.exposure(t,j)     = 0;    
               p.ind_weight(t,j)   = 0;

            end   

        end

    end

end

% Initialize a strucutre str where we store the results at aggregate
% portfolio level
str.caldt        = p.caldt;
str.available    = sum(~isnan(p.ret),2,'omitmissing'); % how many industries were available each day (some industries may became available later in the database)
str.port_weight  = p.ind_weight./str.available; % optimal weight at portfolio level. See Eq.8

% Limit the exposure of each industry at "max_not_trade"
idx_above_max_not = find(str.port_weight > max_not_trade);
str.port_weight(idx_above_max_not) = max_not_trade;

% In case the portfolio exposure exceed the maximum leverage, rescaled each
% position proportionally
str.sum_exposure  = sum(str.port_weight,2,'omitmissing'); % implied daily leverage
idx_above_max_lev = find(str.sum_exposure > max_leverage); % find days when the maximum leverage was exceeded
str.port_weight(idx_above_max_lev,:) = str.port_weight(idx_above_max_lev,:)./str.sum_exposure(idx_above_max_lev)*max_leverage; % Eq10
str.sum_exposure  = sum(str.port_weight,2,'omitmissing');

% Compute portfolio daily returns
str.ret_long     = sum(lag_TS(str.port_weight,1).*p.ret,2,'omitmissing');
str.ret_tbill    = (1-sum(lag_TS(str.port_weight,1),2,'omitmissing')).*p.tbill_ret; % return of cash investment/borrowed

if strcmp(invest_cash,"YES") % include cash returns only when selected
    str.ret_long  = str.ret_long + str.ret_tbill;
end

% Compute the timeseries of the AUM
str.AUM       = AUM_0*cumprod(1 + str.ret_long, 'omitnan');
str.AUM_SPX   = AUM_0*cumprod(1 + p.mkt_ret, 'omitnan'); % Calculating the adjusted AUM for SPX

3.1 Initialize Portfolio Variables and Parameters

The script sets initial conditions including assets under management (AUM), whether cash is invested in short-term bonds, target volatility, maximum leverage, and limits on industry exposure. Arrays for tracking exposure, weights, and trailing stops are initialized based on the number of industries and data points.

Click to see the MATLAB code

MATLAB

AUM_0 = 1;
invest_cash = "YES";
target_vol = 0.02;
max_leverage = 2;
max_not_trade = 0.20;
p.exposure = zeros(size(p.price));
p.ind_weight = zeros(size(p.price));
p.trail_stop_long = NaN(size(p.price));
N_ind = length(p.industry_names);
T = length(p.caldt);

3.2 Calculate Exposure and Weights Dynamically

The nested loops traverse each industry and trading day. The script checks for valid data and calculates exposure based on the presence of a long signal and the status of trailing stops. When conditions are met, exposure and industry weights are adjusted to reflect the desired risk level and strategy rules.

Click to see the MATLAB code

MATLAB

for j = 1:length(p.industry_names)
    for t = 1:T
        if ~isnan(p.ret(t,j)) && ~isnan(p.long_band(t,j))
            if p.exposure(t-1,j) <= 0 && p.long_signal(t,j) == 1
                p.exposure(t,j) = 1;
                p.trail_stop_long(t,j) = p.short_band(t,j);
                lev_vol = target_vol/p.vol(t,j);
                p.ind_weight(t,j) = p.exposure(t,j) * lev_vol;
            elseif p.exposure(t-1,j) == 1 && p.price(t,j) > max(p.trail_stop_long(t-1,j), p.short_band(t,j))
                p.exposure(t,j) = 1;
                p.trail_stop_long(t,j) = max(p.trail_stop_long(t-1,j), p.short_band(t,j));
                lev_vol = target_vol/p.vol(t,j);
                p.ind_weight(t,j) = p.exposure(t,j) * lev_vol;
            elseif p.exposure(t-1,j) == 1 && p.price(t,j) <= max(p.trail_stop_long(t-1,j), p.short_band(t,j))
                p.exposure(t,j) = 0;
                p.ind_weight(t,j) = 0;
            end
        end
    end
}

3.3 Aggregate and Adjust Portfolio Weights

The script adjusts portfolio weights to ensure compliance with risk management rules such as maximum notional trade size per industry and overall leverage limits. Weights are scaled down if the aggregate exposure exceeds predefined thresholds.

Click to see the MATLAB code

MATLAB

str.caldt = p.caldt;
str.available = sum(~isnan(p.ret), 2, 'omitmissing');
str.port_weight = p.ind_weight./str.available;
idx_above_max_not = find(str.port_weight > max_not_trade);
str.port_weight(idx_above_max_not) = max_not_trade;
str.sum_exposure = sum(str.port_weight, 2, 'omitmissing');
idx_above_max_lev = find(str.sum_exposure > max_leverage);
str.port_weight(idx_above_max_lev, :) = str.port_weight(idx_above_max_lev, :) ./ str.sum_exposure(idx_above_max_lev) * max_leverage;

3.4 Compute Portfolio Returns

Calculates daily returns for the portfolio, including returns from cash investments if specified. The returns are then used to compute the growth of assets under management (AUM), both for the strategy and relative to a market index, providing a measure of performance over time.

Click to see the MATLAB code

MATLAB

str.ret_long = sum(lag_TS(str.port_weight, 1) .* p.ret, 2, 'omitmissing');
str.ret_tbill = (1 - sum(lag_TS(str.port_weight, 1), 2, 'omitmissing')) * p.tbill_ret;
if strcmp(invest_cash, "YES")
    str.ret_long = str.ret_long + str.ret_tbill;
endif
str.AUM = AUM_0 * cumprod(1 + str.ret_long, 'omitnan');
str.AUM_SPX = AUM_0 * cumprod(1 + p.mkt_ret, 'omitnan');

4. Study Results

Overview

Step 4 involves plotting and showing the stats of the backtest.

Step 4 Code

Click to see the MATLAB code

MATLAB

freq = 21;
fig = figure(); 
plot(str.caldt([1:freq:end end]), str.AUM([1:freq:end end]), 'LineWidth', 2,'Color','k'); hold on
plot(str.caldt([1:freq:end end]), str.AUM_SPX([1:freq:end end]), 'LineWidth', 1, 'Color','r'); hold off
grid off; set(gca,'GridLineStyle',':');
set(gca, 'GridLineStyle', ':');
set(gca, 'XTick');
xlim([str.caldt(1) str.caldt(end)]);
datetick('x', 'yyyy', 'keepticks'); % Keep the ticks and format them

xtickangle(90); ytickformat('$%,.0f');
ax=gca; ax.YAxis.Exponent = 0;
set(gca, 'FontSize', 8); % Adjust font size
legend('Timing Ind.','Market','Location','northwest');
title(['Timing ' num2str(N_ind) ' Industries Strategy'], 'FontWeight', 'bold','FontSize',12);

y_tick(1) = AUM_0;
for j = 2:1:100, y_tick(j) = y_tick(j-1)*4; end
yticks([y_tick]); 

set(gca, 'YScale', 'log');

table_statistics = table_of_stats([str.AUM str.AUM_SPX],[size(p.ret,2) 1],{"Timing Ind.","Market"},str.caldt,p.mkt_ret,p.tbill_ret);
display(table_statistics)

Needed functions for Step 4:

table_of_stats
SortinoRatio
computeHitRatio
computeSkewness

Click to see the MATLAB code of table_of_stats

MATLAB

function table_statistics = table_of_stats(AUM,N_assets,names,caldt,mkt_ret,tbill_ret)

    table_statistics = struct();

    for s = 1:size(AUM,2)
        
        ret = price2return(AUM(:,s),1);

        table_statistics(s).strategy = names{s};
        
        table_statistics(s).assets   = N_assets(s);
        table_statistics(s).irr      = round((prod(1+ret,'omitmissing')^(252/length(ret))-1)*100,1);
        table_statistics(s).vol      = round(std(ret,'omitmissing')*sqrt(252)*100,1);
        table_statistics(s).sr       = round(mean(ret,'omitmissing')/std(ret,'omitmissing')*sqrt(252),2);
        table_statistics(s).sortino  = round(SortinoRatio(ret)*sqrt(252),2);
        
        table_statistics(s).ir_d     = computeHitRatio(ret,caldt,"daily");
        table_statistics(s).ir_m     = computeHitRatio(ret,caldt,"monthly");
        table_statistics(s).ir_y     = computeHitRatio(ret,caldt,"yearly");

        table_statistics(s).skew_d   = round(skewness(ret),2);
        table_statistics(s).skew_m   = round(computeSkewness(ret,caldt,"monthly"),2);
        table_statistics(s).skew_y   = round(computeSkewness(ret,caldt,"yearly"),2);
        table_statistics(s).mdd      = round(maxdrawdown(AUM(:,s))*100,0);

        XX            = [mkt_ret-tbill_ret];
        YY            = ret-tbill_ret;
        
        regr         = fitlm(XX,YY);
        coef         = regr.Coefficients.Estimate;
        tstas        = regr.Coefficients.tStat;

        table_statistics(s).alpha    = round(coef(1)*100*252,2);
        table_statistics(s).alpha_t  = tstas(1);

        table_statistics(s).beta     = round(coef(2),2);
        table_statistics(s).beta_t   = tstas(2);

        table_statistics(s).worst_ret= round(min(ret)*100,2);
        table_statistics(s).worst_day= datestr(caldt(find(ret == min(ret))));
        table_statistics(s).best_ret = round(max(ret)*100,2);
        table_statistics(s).best_day = datestr(caldt(find(ret == max(ret))));
                
    end
    
    originalTable = struct2table(table_statistics);

    % Step 1: Extract the column names
    colNames = originalTable.Properties.VariableNames;

    % Step 2: Transpose the data
    transposedData = table2array(originalTable)';

    % Step 3: Create the new table
    newTable = array2table(transposedData);%, 'VariableNames', colNames);

    % Step 4: Add the original column names as a new column
    newTable = addvars(newTable, colNames', 'Before', 1);
    
    % Step 5: Name the column of the new Table
    newTable.Properties.VariableNames = [" " names];
    

    % Display the result
    table_statistics = newTable;
    
end

Click to see the MATLAB code of SortinoRatio

MATLAB

function sortinoRatio = SortinoRatio(returns)
    % Compute the mean return
    meanReturn = mean(returns,'omitmissing');

    % Compute the downside deviation
    downsideReturns   = returns(returns < 0);
    downsideDeviation = std(downsideReturns,'omitmissing');

    % Compute the Sortino Ratio
    sortinoRatio = (meanReturn) / downsideDeviation;
end

Click to see the MATLAB code of computeHitRatio

MATLAB

function hitRatio = computeHitRatio(ret,caldt, timeframe)
    % Ensure the timeframe is valid
    validTimeframes = {'daily', 'monthly', 'yearly'};
    if ~ismember(timeframe, validTimeframes)
        error('Invalid timeframe. Choose "daily", "monthly", or "yearly".');
    end

    % Convert datenum to datetime for easier manipulation
    dates   = datetime(caldt, 'ConvertFrom', 'datenum');
    returns = ret;

    % Calculate the Hit Ratio based on the timeframe
    switch timeframe
        case 'daily'
            hitRatio = round(sum(returns > 0) / sum(abs(returns) > 0) * 100, 0);

        case 'monthly'
            % Group by year and month
            [yearMonth, ~, idx] = unique(dates.Year*100 + dates.Month);
            monthlyReturns = accumarray(idx, returns, [], @(x) prod(1 + x,'omitmissing') - 1);
            hitRatio = round(sum(monthlyReturns > 0) / sum(abs(monthlyReturns) > 0) * 100, 0);

        case 'yearly'
            % Group by year
            [years, ~, idx] = unique(dates.Year);
            yearlyReturns = accumarray(idx, returns, [], @(x) prod(1 + x,'omitmissing') - 1);
            hitRatio = round(sum(yearlyReturns > 0) / sum(abs(yearlyReturns) > 0) * 100, 0);

        otherwise
            error('Unexpected timeframe.');
    end
end

Click to see the MATLAB code of computeSkewness

MATLAB

function skw = computeSkewness(ret,caldt, timeframe)
    % Ensure the timeframe is valid
    validTimeframes = {'daily', 'monthly', 'yearly'};
    if ~ismember(timeframe, validTimeframes)
        error('Invalid timeframe. Choose "daily", "monthly", or "yearly".');
    end

    % Convert datenum to datetime for easier manipulation
    dates = datetime(caldt, 'ConvertFrom', 'datenum');
    returns = ret;

    % Calculate the Hit Ratio based on the timeframe
    switch timeframe
        case 'daily'
            skw = skewness(returns);

        case 'monthly'
            % Group by year and month
            [yearMonth, ~, idx] = unique(dates.Year*100 + dates.Month);
            monthlyReturns = accumarray(idx, returns, [], @(x) prod(1 + x,'omitmissing') - 1);
            skw = skewness(monthlyReturns);

        case 'yearly'
            % Group by year
            [years, ~, idx] = unique(dates.Year);
            yearlyReturns = accumarray(idx, returns, [], @(x) prod(1 + x,'omitmissing') - 1);
            skw = skewness(yearlyReturns);;

        otherwise
            error('Unexpected timeframe.');
    end
end

4.1 Compute Portfolio Returns

This script visualizes the performance of the investment strategy relative to a market index. It plots the cumulative asset value of the strategy and the market index on a logarithmic scale to better illustrate growth trends over time. The graph is customized with specific formatting options to enhance readability and presentation.

Click to see the MATLAB code

MATLAB

freq = 21;
fig = figure();
plot(str.caldt([1:freq:end end]), str.AUM([1:freq:end end]), 'LineWidth', 2, 'Color', 'k'); hold on
plot(str.caldt([1:freq:end end]), str.AUM_SPX([1:freq:end end]), 'LineWidth', 1, 'Color', 'r'); hold off
grid off; set(gca, 'GridLineStyle', ':');
set(gca, 'XTick');
xlim([str.caldt(1) str.caldt(end)]);
datetick('x', 'yyyy', 'keepticks');
xtickangle(90); ytickformat('$%,.0f');
ax = gca; ax.YAxis.Exponent = 0;
set(gca, 'FontSize', 8);
legend('Timing Ind.', 'Market', 'Location', 'northwest');
title(['Timing ' num2str(N_ind) ' Industries Strategy'], 'FontWeight', 'bold', 'FontSize', 12);
y_tick(1) = AUM_0;
for j = 2:1:100, y_tick(j) = y_tick(j-1)*4; end
yticks([y_tick]);
set(gca, 'YScale', 'log');

4.2 Compute and Display Summary Statistics

Generates a summary statistics table that includes important metrics such as internal rates of return, volatility, Sharpe ratio, and other performance indicators. This table is calculated using a custom function table_of_stats, which processes the asset management data and returns to provide a comprehensive statistical breakdown of the strategy’s performance over time.

Click to see the MATLAB code

MATLAB

table_statistics = table_of_stats([str.AUM str.AUM_SPX], [size(p.ret,2) 1], {"Timing Ind.", "Market"}, str.caldt, p.mkt_ret, p.tbill_ret);
display(table_statistics)

Full Code and Parameters Editing

You can download the full code from here and run it directly cell by cell in your MATLAB!
Remember to run cell 0 (at the bottom) to load needed functions!

Click on the following Link to download the file:

Download Backtesting_A-Century-of-Industry-Trends

If you are interested in a version where you only want to change parameters, this is a version where the backtest is put in a function and you just edit the parameters to see the results faster.

Download Backtesting_A-Century-of-Industry-Trends (Easy Parameter)

If you have questions feel free to contact us on twitter/X or contact us directly here

Teaming up with well-known researcher and industry expert @GaryAntonacci , we are proud to announce the publication of our new research piece

A Century of Profitable Industry Trends

This study delves into the profitability of a simple, low-frequency trend-following model… pic.twitter.com/1j7etbXdHq
— Concretum Research (@ConcretumR) June 7, 2024