function [OptimalModelParameters, optimalmodelorder, ModelOutput, A, R2] = optimalARmodelfit(y)

% Inputs:
% y: observed data (dependent variable); it must be a column vector
% Outputs:
% OptimalModelParameters: vector of estimated optimal parameter values
% optimalmodelorder: optima AR model order estimated by AIC
% y_fit: fitted y's obtained as output of the optimal model
% A: optimal regressor matrix
% R2: coefficient of determination of the optimal model

y = y(:);

for n = 1:40
A = zeros(length(y)-n,n);
for i = n+1:length(y)
    A(i-n,:) = y(i-1:-1:i-n)';
end
b = y(n+1:end);
theta = A\b;
y_fit = A*theta;
AIC(n) = length(y)*log(1/length(y)*sum((y(n+1:end)-y_fit).^2))+2*(n);
if n>2
    if AIC(n-2)>AIC(n-1) && AIC(n)>AIC(n-1)
        break
    end
end
end
optimalmodelorder = n-1;

clear A y_fit
for i = optimalmodelorder+1:length(y)
    A(i-optimalmodelorder,:) = y(i-1:-1:i-optimalmodelorder)';
end
b = y(optimalmodelorder+1:end);
OptimalModelParameters = A\b;
ModelOutput = A*OptimalModelParameters;
% figure,
% plot((1:length(y)-optimalmodelorder),y(optimalmodelorder+1:end),'.r',(1:length(y)-optimalmodelorder),ModelOutput,'-k','LineWidth',1,'MarkerSize',14)
% xlabel('x')
% ylabel('y')
% title( sprintf( 'Observations and predicted outcome from an autoregressive model of order %f', optimalmodelorder ) );
% legend('Data','Fitted model')

R2 = 100*(1-sum((y(optimalmodelorder+1:end)-ModelOutput).^2)/sum((y(optimalmodelorder+1:end)-mean(y(optimalmodelorder+1:end))).^2));

