Needs["MultivariateStatistics`"]

Stations = {"Munich","Cape Town"};

data=WeatherData[#,"MeanTemperature",{{2005,1,1},{2009,1,1},"Day"}]&/@ Stations;

(*Only dates that are in both lists*)
dates=Intersection[data[[1,All,1]],data[[2,All,1]]];

PointsETSE=Cases[data[[1]],{Alternatives@@dates,_}]; (* data from first item in Stations*)

PointsETSL = Cases[data[[2]],{Alternatives@@dates,_}]; (*data from second item in Stations*)
T1 = PointsETSE[[All,2]];
T2 =PointsETSL[[All,2]];
Pairs=Table[{T1,T2}[[All,i]],{i,1,Length[T1]}];
GaltonDat = Pairs;
GaltonX =T1;
GaltonY = T2;

Scatterplot

legend=Framed@Column[
Table[
Style[
Stations[[i]],ColorData[3,i],Bold,Background->White,"Graphics",15
],{i,Length[Stations]}
]
];

DateListPlot[data
,Epilog->Inset[legend,{{2008,7,21},-5}]
,PlotStyle->ColorData[3,"ColorList"]
,BaseStyle->{FontFamily->"Times",FontSize->14}
,FrameLabel->{"time","temperature [°C]"}
]
GaltonScatter=ListPlot[Pairs
,AspectRatio->35/50
,BaseStyle->{FontFamily->"Times",FontSize->14}
,Joined->False
,PlotRange->{{-15,35},{0,35}}
,Frame->True
,FrameLabel->{"temperature in Munich [°C]","temperature in Cape Town [°C]"}
,FrameTicks->{LinTicks[-20,30,5,2],LinTicks[-20,30,5,2],LinTicks[-20,30,5,2,ShowTickLabels->False],LinTicks[-20,30,5,2,ShowTickLabels->False]}
]

some descriptive Statistics

{GalX,GalY}=Transpose[GaltonDat];
{{minx,maxx},{miny,maxy}}=Map[{Min[#],Max[#]}&,{GalX,GalY}]

SDevX=StandardDeviation[GaltonX]
SDevY=StandardDeviation[GaltonY]
MY=Mean[GaltonY]
MX=Mean[GaltonX]

Fit  regression

some variations ... all result in the same thing

(*lmTr=LinearModelFit[Transpose@{GaltonX,GaltonY},{1,x},x];
lmTr["BestFit"]*)
LinRegr=Fit[GaltonDat,{1,x},x]
lm=LinearModelFit[GaltonDat,{1,x},x]
lm["BestFit"]

q=EllipsoidQuartiles[GaltonDat];

"LinearModelFit" offers some more options :

lm[{"RSquared","ANOVATable"}]

some bivariate statistics

r2=Correlation[GaltonX,GaltonY]

CovXY=Covariance[GaltonX,GaltonY]
CovYX=Covariance[GaltonY,GaltonX]

r=(CovXY/(SDevX*SDevY))

constructing the regression line manually

slope1=r2*SDevY/SDevX
intercept1=MY-(slope1*MX)
regrLine=intercept1+slope1*x

constructing the SDLine manually

slope2=Sign[r2]*SDevY/SDevX
intercept2=(MY-(slope2*MX))
SDLine=slope2*x+intercept2

Discussion of plot

regrLine is the same as LinRegr (green and blue)
SD line is red

OneSDrightofMean=MX+SDevX
OneSDleftofMean = MX-SDevX

LinesCPH=Plot[{regrLine,SDLine,LinRegr},{x,-20,40}
,AspectRatio->1
,PlotRange->{{0,80},{0,80}}
,PlotStyle->{
{RGBColor[0,1,0],Thickness[0.01]} (*green*)
,{RGBColor[1,0,0],Thickness[0.01]} (*red*)
,{RGBColor[0,0,1],Thickness[0.0025]}(*blue*)
}
];

PlotOneSDRightofFathersMean=ListLinePlot[{{OneSDleftofMean,0},{OneSDleftofMean,40}}
,PlotStyle->{Orange,Dashed,Thick}];
PlotOneSDTopOfSonsMean=Plot[MY+SDevY,{x,-20,40}
,PlotStyle->{Yellow,Dashed,Thick}];
Show[GaltonScatter,
Graphics[{GrayLevel[0.8],Thickness[0.01],q[[1]]}],
Graphics[{GrayLevel[0.5],Thickness[0.005],q[[2]]}],
Graphics[{GrayLevel[0.1],Thickness[0.001],q[[3]]}],
LinesCPH,
PlotOneSDRightofFathersMean
,PlotOneSDTopOfSonsMean]