I’m working in MATLAB to create a value stock portfolio using historical series through the Fama-French 5-factor model. As you can see, I’ve used historical data from 1999 to 2005 to identify (from a random sample of S&P 500 stocks) value stocks, which I then used to create a portfolio using Markowitz optimization, with the standard constraints of non-negative weights and the sum of weights equal to 1.

The problem arises during the portfolio rotation: my idea is to re-run the FF5 analysis every two years to identify new value stocks to add to the portfolio, and so on for each year of analysis. Where am I going wrong? Thanks in advance to everyone.

Here is my code:

% Caricamento dei dati da Excel
priceData = readtable('PRICES-1.xlsx', 'VariableNamingRule', 'preserve');
ff5Data = readtable('FF5.xlsx', 'VariableNamingRule', 'preserve');
forecastData = readtable('FORECAST.xlsx', 'VariableNamingRule', 'preserve');
ff5ForecastData = readtable('FF5_FORECAST.xlsx', 'VariableNamingRule', 'preserve');

% Estrarre le date e i rendimenti dei titoli da priceData
dates = priceData{:, 1};
returns = priceData{:, 2:end};
titleNames = priceData.Properties.VariableNames(2:end);

% Convertire returns in numerico se necessario
if iscell(returns)
    returns = cellfun(@str2double, returns);
end

% Estrarre i fattori di Fama e French e il tasso privo di rischio
rmrf = ff5Data{:, 'Mkt-RF'};
smb = ff5Data{:, 'SMB'};
hml = ff5Data{:, 'HML'};
rmw = ff5Data{:, 'RMW'};
cma = ff5Data{:, 'CMA'};
rf = ff5Data{:, 'RF'};

% Convertire le colonne dei fattori in numerico (se necessario)
factor_columns = {rmrf, smb, hml, rmw, cma, rf};
for i = 1:length(factor_columns)
    if iscell(factor_columns{i})
        factor_columns{i} = cellfun(@str2double, factor_columns{i});
    end
end
[rmrf, smb, hml, rmw, cma, rf] = deal(factor_columns{:});

% Numero di titoli
num_titoli = size(returns, 2);

% Matrice dei fattori senza l'intercetta
factors = [rmrf, smb, hml, rmw, cma];

% Matrice per i coefficienti di regressione
coefficients = zeros(num_titoli, 6);

% Regressione robusta per ogni titolo
for i = 1:num_titoli
    Y = returns(:, i) - rf; % Rendimento in eccesso rispetto al tasso privo di rischio
    validIdx = ~isnan(Y) & all(~isnan(factors), 2); % Indici validi (senza NaN)
    B = robustfit(factors(validIdx, :), Y(validIdx));
    coefficients(i, :) = B'; % I coefficienti includono l'intercetta come primo elemento
end

% Identifica i titoli "value" usando il coefficiente HML
threshold = prctile(coefficients(:, 5), 70); % Calcola il 70° percentile dei coefficienti HML
value_titles = coefficients(:, 5) > threshold;

% Visualizza i risultati con i nomi dei titoli
value_titles_names = titleNames(value_titles);
disp('Titoli "value":');
disp(strjoin(value_titles_names, ', '));

% Caricare i dati da FORECAST.xlsx
forecastPrices = forecastData{:, 2:end};
forecastDates = forecastData{:, 1};
forecastTitleNames = forecastData.Properties.VariableNames(2:end);

% Convertire i prezzi in numerico (se necessario)
if iscell(forecastPrices)
    forecastPrices = cellfun(@str2double, forecastPrices);
end

% Filtrare i titoli "value" nei dati di forecast
[~, valueIdx] = ismember(value_titles_names, forecastTitleNames);
valueForecastPrices = forecastPrices(:, valueIdx);

% Calcolare i rendimenti logaritmici dai dati di forecast
valueReturns = log(valueForecastPrices(2:end, :) ./ valueForecastPrices(1:end-1, :));
forecastDates = forecastDates(2:end);

% Costruzione del portafoglio con il modello di Markowitz
num_value_titles = size(valueReturns, 2);
mu = mean(valueReturns)'; % Rendimenti attesi
Sigma = cov(valueReturns); % Matrice di covarianza

% Ottimizzazione di Markowitz con vincoli
Aeq = ones(1, num_value_titles);
beq = 1;
lb = zeros(num_value_titles, 1);

% Rolling window parameters
window_size = 252; % Dimensione della finestra rolling (giornaliera, circa 1 anno)
num_periods = size(valueReturns, 1) - window_size;
rolling_returns = zeros(num_periods, 1);

options = optimoptions('quadprog', 'Display', 'off');

% Periodo iniziale di due anni
initial_period = 2 * 252; % Assumendo 252 giorni di trading all'anno sui 365

% Calcolare i rendimenti per i primi due anni
for t = 1:initial_period
    window_returns = valueReturns(t:t+window_size-1, :);
    mu_rolling = mean(window_returns)';
    Sigma_rolling = cov(window_returns);
    [w, ~] = quadprog(Sigma_rolling, [], [], [], Aeq, beq, lb, [], [], options);
    rolling_returns(t) = w' * valueReturns(t+window_size, :)';
end

% Iterare ogni due anni
current_period = initial_period + 1;
while current_period <= num_periods
    % Dati di prezzo e fattori per i due anni precedenti
    past_returns = valueReturns(current_period-initial_period:current_period-1, :);
    past_factors = ff5ForecastData{current_period-initial_period:current_period-1, 2:end};

    % Convertire i fattori in numerico se necessario
    if iscell(past_factors)
        past_factors = cellfun(@str2double, past_factors);
    end

    % Regressione robusta per ogni titolo
    past_coefficients = zeros(num_titoli, 6);
    for i = 1:num_titoli
        Y = past_returns(:, i) - rf(current_period-initial_period:current_period-1);
        validIdx = ~isnan(Y) & all(~isnan(past_factors), 2); % Indici validi (senza NaN)
        B = robustfit(past_factors(validIdx, :), Y(validIdx));
        past_coefficients(i, :) = B';
    end
    
    % Identifica nuovi titoli "value"
    threshold = prctile(past_coefficients(:, 5), 70);
    new_value_titles = past_coefficients(:, 5) > threshold;
    new_value_titles_names = titleNames(new_value_titles);
    
    % Aggiungi i nuovi titoli "value" al portafoglio
    [~, new_valueIdx] = ismember(new_value_titles_names, forecastTitleNames);
    valueIdx = unique([valueIdx; new_valueIdx]);
    valueForecastPrices = forecastPrices(:, valueIdx);
    valueReturns = log(valueForecastPrices(2:end, :) ./ valueForecastPrices(1:end-1, :));
    
    % Calcolare i rendimenti logaritmici dai dati di forecast
    num_value_titles = size(valueReturns, 2);
    mu = mean(valueReturns)'; % Rendimenti attesi
    Sigma = cov(valueReturns); % Matrice di covarianza
    
    % Calcolare i rendimenti per i successivi due anni
    for t = current_period:current_period + 2 * 252 - 1
        if t + window_size > size(valueReturns, 1)
            break;
        end
        window_returns = valueReturns(t:t+window_size-1, :);
        mu_rolling = mean(window_returns)';
        Sigma_rolling = cov(window_returns);
        [w, ~] = quadprog(Sigma_rolling, [], [], [], Aeq, beq, lb, [], [], options);
        rolling_returns(t) = w' * valueReturns(t+window_size, :)';
    end
    
    % Avanza al periodo successivo
    current_period = current_period + 2 * 252;
end

% Converti i rendimenti rolling in percentuale
rolling_returns_percentage = rolling_returns * 100;

% Output dei rendimenti rolling del portafoglio
disp('Rendimenti rolling del portafoglio (%):');
disp(rolling_returns_percentage);

% Calcolo della performance complessiva del portafoglio
cumulative_returns = cumprod(1 + rolling_returns) - 1;
cumulative_returns_percentage = cumulative_returns * 100;

disp('Rendimento cumulativo del portafoglio (%):');
disp(cumulative_returns_percentage(end));

The first part is ok, is the iterative process that does not work!