Practical Linear Algebra: A Geometry Toolbox, Third edition
Gerald Farin & Dianne Hansford
CRC Press, Taylor & Francis Group, An A K Peters Book, 2013

Least Squares Approximation / Linear Regression
Find the least squares best fit line to data.
Use data from Mathematica’s Financial Data and use the best fit line for spotting a trend in the data. 
This example is for AIG’s stock from 1/1/2000 to today.
Also: input test data sets
This nb uses the normal equations and Mathematica’s LinearSolve routine.
See PLA Section 12.7 Gauss for Linear Systems, Least Squares

Initialize output folders and graphics

(* to enable file export: *)
(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];)
SetOptions[Plot,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];
SetOptions[ListPlot,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];
SetOptions[DateListPlot,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];

Create data

list=FinancialData["AP","Jan. 1, 2000"];
(* Also: GM, AIG, AP *)

(* Print["financial data = ",MatrixForm[list]]; *)

length=Length[list];
data=Table[list[[i,2]],{i,1,length}];
datplot=ListPlot[data,PlotStyle->{GrayLevel[0.4],PointSize[0.004]}];

mat=Table[{i,1},{i,1,length}];
(* Print["mat=",MatrixForm[mat]]; *)

Solve least squares and form line for graphics

data:  (Subscript[x, i]), (Subscript[y, i]),  i=1,...,N

want
	y = ax + b

such that 

	(Subscript[y, 1]) = a (Subscript[x, 1]) + b
	⋮
	⋮
	(Subscript[y, N]) = a (Subscript[x, N]) + b
	
yields overdetermined linear system for a = a,b

	(({{(Subscript[y, 1])}, {⋮}, {⋮}, {(Subscript[y, N])}})) =  a (({{(Subscript[x, 1])}, {⋮}, {⋮}, {(Subscript[x, N])}}))  + (({{b}, {⋮}, {⋮}, {b}})).    ⟹  (({{(Subscript[y, 1])}, {⋮}, {⋮}, {(Subscript[y, N])}})) = (({{(Subscript[x, 1]), 1}, {({{⋮}, {⋮}}), ({{⋮}, {⋮}})}, {(Subscript[x, N]), 1}}))(({{a}, {b}}))
	
Short:

	Y = Xa   	where (X^T)X is a 2x2 matrix. These are the normal equations.
	
	(X^T)Y = (X^T)X a	
	
Solution
	a = ((((X^T)X))^(-1))Y

(*Solve the least squares system: *)
Print["Least squares line coefficients:"];
coeffs=LinearSolve[Transpose[mat].mat,Transpose[mat].data]

(*the least squares line line[x]: *)
line[x_]:=coeffs[[1]]*x+coeffs[[2]];
lineplot=Plot[line[x],{x,1,length},PlotStyle->{Thickness[0.01], Black}];

(* pic=Show[{datplot,lineplot}] *)


(* ---------Better than above plot: uses date in plot --------- *)
(* the data run from 1 ... length; since they are dates, use the integer i as the x-coordinate *)
linelist = {{list[[1,1]],line[1]},{list[[length,1]],line[length]}} 

pic =DateListPlot[{list,linelist},Joined->{False,True}, GridLines->False,Frame->True,
FrameTicks->{{Automatic,None},{Automatic,None}},PlotStyle->{{GrayLevel[0.5],PointSize[0.005]},{GrayLevel[0], Thickness[0.009]}}]
Export["lsq_Dow.eps",pic];

More data sets -- repeat least squares best fit line

(* For example in gensys *)
experiment1 = {{0,30},{10,25},{20,40},{30,40},{40,30},{50,5},{60,25}};
(* For exercises *)
experiment2 = {{-2,2},{-1,1},{0,0},{1,1},{2,2}};
experiment3 = {{0,0},{1,1},{2,2},{1,-1},{2,-2}};
experiment4 = {{-4,1},{4,3}}

experiment = experiment2;

elngth = Length[experiment];
emat = Table[{experiment[[i,1]],1},{i,1,elngth}];
edata = Table[experiment[[i,2]],{i,1,elngth}];

Print["matrix =",MatrixForm[emat],"  rhs=",MatrixForm[edata]];

(*Solve the least squares system: *)
ecoeffs=LinearSolve[Transpose[emat].emat,Transpose[emat].edata];

Print["normal equations: ",MatrixForm[Transpose[emat].emat],"u =",MatrixForm[Transpose[emat].edata]];
Print["coefficients =",ecoeffs,"  as decimal ",ecoeffs//N];

eline[x_]:=ecoeffs[[1]]*x+ecoeffs[[2]];

elinelist = {{experiment[[1,1]],eline[experiment[[1,1]]]},
{experiment[[elngth,1]],eline[experiment[[elngth,1]]]}};
Print["Line approx: end1=",elinelist[[1]],"  end2=",elinelist[[2]]];


pic2 = ListPlot[{experiment,elinelist},Joined->{False,True},PlotStyle->{{GrayLevel[0.5],PointSize[0.02]},{GrayLevel[0],Thickness[0.008]}}]
Export["lsq_experiment.eps",pic2];