function Treg = make_regression_table(tbl, delayVar, p, horizon)
% Build regression table:
% Inputs: delayLag1...delayLagP + current exogenous features
% Output: targetDelay = delay(t+horizon)

if nargin < 4
    horizon = 1;
end

delay = tbl.(delayVar);
N = height(tbl);

workTbl = tbl;

for k = 1:p
    lagName = sprintf('delayLag%d', k);
    workTbl.(lagName) = [nan(k,1); delay(1:end-k)];
end

workTbl.targetDelay = [delay(1+horizon:end); nan(horizon,1)];

validRows = true(N,1);
for k = 1:p
    lagName = sprintf('delayLag%d', k);
    validRows = validRows & ~isnan(workTbl.(lagName));
end
validRows = validRows & ~isnan(workTbl.targetDelay);

allNames = workTbl.Properties.VariableNames;
removeNames = [{delayVar}, {'targetDelay'}];
featureNames = allNames(~ismember(allNames, removeNames));

orderedLagNames = cell(1,p);
for k = 1:p
    orderedLagNames{k} = sprintf('delayLag%d', k);
end

exoNames = featureNames(~startsWith(featureNames, 'delayLag'));
finalOrder = [orderedLagNames, exoNames, {'targetDelay'}];

Treg = workTbl(validRows, finalOrder);
end
