Hi Anton

2) I don't agree , the cloud of points plot is a Poincaré section of some periodic dynamic signal, and the symmetry is not mandatory.

For example consider :

ListPlot[Partition[Table[Sin[t]/t, {t, 1, 100, 1}], 2, 1]]

which gives a kind of multispiral plot.

3) I did some trials will the Clayton copula, not so sure now it would give the best fit. But quite sure that the copula approach is worth to try.

The Clayton copula has higher correlation for the head of marginals but Mathematica does not have the inverted Clayton availaible. The Gumbel correlates both head and tail ends, and it is easy to use it in Mathematica.

Here is the idea with the Gumbel :

(I am not used to do fitting of data, so this just gives an idea and this is far to be a good fit)

dist = CopulaDistribution[{"GumbelHougaard", 4}, {NormalDistribution[16, 8], NormalDistribution[16, 8]}];

ListPlot[RandomVariate[dist, 3000], PlotRange -> {{-10, 35}, {-10, 35}}]

which gives for 0° a plot which shows that the most probable temperature is higher, this is coherent with your results

Plot[Evaluate[PDF[dist, {Ty, Tt}]]/PDF[NormalDistribution[16, 8], Ty] /. Ty -> 0, {Tt, -5, 10}]

The overall behaviour of the distribution : symmetry in the middle range, and skew trend in the high temperatures seems to be correct,

nevertheless the model is very approximative and far to be as much accurate than with your quantile regression package.