% Regression

%% 
% Generating the model
% y = X*B + eps;
clear;
clc;
x = load('regress_data.txt');
y = x(:,1);
n = length(y);
X = [ones(n,1) x(:,2:end-1)];


%%
% Using the command regress from the statistics toolbox we estimate the
% betas's as follow

[b_st,bint,r,rint,stats] = regress(y,X);
sighat = sqrt(stats(4));

%%
% using the command fminunc from the optimization toolbox we estimate the
% betas's as follow

opt = optimset('Display','iter','MaxFunEvals',400,'MaxIter',400,'TolFun',1.0000e-010,...
    'TolX',1.0000e-006,'FunValCheck','off','OutputFcn',[],'ActiveConstrTol',[],'BranchStrategy',...
    [],'DerivativeCheck','off','Diagnostics','off','DiffMaxChange',0.1000,'DiffMinChange',1.0000e-008,...
    'GoalsExactAchieve',[],'GradConstr',[],'GradObj','on','Hessian','off','HessMult',[],'HessPattern',...
    'sparse(ones(numberofvariables))','HessUpdate','bfgs','InitialHessType','scaled-identity','InitialHessMatrix',...
    [],'Jacobian',[],'JacobMult',[],'JacobPattern',[],'LargeScale','on','LevenbergMarquardt',[],'LineSearchType',...
    [],'MaxNodes',[],'MaxPCGIter','max(1,floor(numberofvariables/2))','MaxRLPIter',[],'MaxSQPIter',[],'MaxTime',[],...
    'MeritFunction',[],'MinAbsMax',[],'NodeDisplayInterval',[],'NodeSearchStrategy',[],'NonlEqnAlgorithm',[],...
    'NoStopIfFlatInfeas',[],'PhaseOneTotalScaling',[],'Preconditioner',[],'PrecondBandWidth',0,'RelLineSrchBnd',[],...
    'RelLineSrchBndDuration',[],'ShowStatusWindow',[],'Simplex',[],'TolCon',[],'TolPCG',0.1000,'TolRLPFun',[],...
    'TolXInteger',[],'TypicalX','ones(numberofvariables,1)');
warning off all;
Ini_Guess = [0.1 0.3 0.2 0.5 0.9 1];
[par,fval,exitflag,output] = fminunc('OLS_b',Ini_Guess,opt,X,y)

% disp('beta St. Tool');
% disp([b_st' sighat]);
% disp('beta Op. Tool');
% disp(par);
% disp(output);

% plot1 = plot(X(:,2),y);
% hold on
% set(gcf,'Color',[1,1,1]);
% set(plot1,'Marker','.','LineStyle','none','Color','r');
% title('Simple Linear Regression','FontWeight','bold');
% xlabel('X','FontWeight','bold');
% ylabel('Y','FontWeight','bold','Rotatio',0.0);
% hold off


