function bfCoefficients = FindPolicyLS(sPaths,T,r,payoffFun,bfHandles)
[numScenarios, numSteps, ~] = size(sPaths);
dt = T/numSteps;
gamma = exp(-r*dt);
numBasis = length(bfHandles); % number of basis functions
bfCoefficients = zeros(numBasis,numSteps-1);
regrMat = zeros(numScenarios, numBasis);
nextValues = zeros(numScenarios, 1);
intrinsicValues = zeros(numScenarios, 1);
% The value function for the last step is exactly the option payoff
% The function squeeze is used to transform a three-dimensional array
% into a one-dimensional vector when we fix two indexes
for s = 1:numScenarios
    nextValues(s) = feval(payoffFun, squeeze(sPaths(s,numSteps,:)));
end
for t = numSteps-1:-1:1
    for s = 1:numScenarios
        prices = squeeze(sPaths(s,t,:));
        intrinsicValues(s) = feval(payoffFun, prices);
        for k = 1:numBasis
            regrMat(s,k) = feval(bfHandles{k}, prices);
        end
    end
    bfCoefficients(:,t) = regrMat\(gamma*nextValues);
    continuationValues = regrMat*bfCoefficients(:,t);
    nextValues = max(continuationValues, intrinsicValues);
end 
