SetDirectory[NotebookDirectory[]]
rawdata=Import["data_for_linear_fit.dat"];
n=Length[rawdata];
c11=n;
M={{c11,c12},{c21,c22}};

"Calculate the best fit parameters using least square method in Matrix form for a set of linear data"
SetDirectory[NotebookDirectory[]]
rawdata=Import["data_for_linear_fit.dat"];
lengthrawdata=Length[rawdata];
xdat[i_]:=rawdata[[i,1]];ydat[i_]:=rawdata[[i,2]];
sumx=Sum[xdat[i],{i,1,lengthrawdata}];
sumy=Sum[ydat[i],{i,1,lengthrawdata}];
sumxsq=Sum[xdat[i]^2,{i,1,lengthrawdata}];
sumysq=Sum[ydat[i]^2,{i,1,lengthrawdata}];
sumxy=Sum[xdat[i]*ydat[i],{i,1,lengthrawdata}];

M={{lengthrawdata,sumx},{sumx,sumxsq}};
X={sumy,sumxy};
{a,b}=Inverse[M].X
lp=ListPlot[rawdata];
plt=Plot[b*x+a,{x,-20,20}];
Show[{lp,plt}]

miny=sumy/lengthrawdata;minx=sumx/lengthrawdata;
ssyy=Sum[(ydat[i]-miny)^2,{i,1,lengthrawdata}];
ssxx=Sum[(xdat[i]-minx)^2,{i,1,lengthrawdata}];
ssxy=Sum[(xdat[i]-minx)(ydat[i]-miny),{i,1,lengthrawdata}];
s=Sqrt[(ssyy-(ssxy^2/ssxx))/(lengthrawdata-1)];
sea=s*Sqrt[(1/lengthrawdata)+(minx^2/ssxx)]
seb=s/Sqrt[ssxx]

"Checking the correctness of standard error in a and b using ParameterTable"
lmf=LinearModelFit[rawdata,{x},x]
lmf[{"BestFit","ParameterTable"}]