Using GLS to fix heteroscedasticity

Question

I have a dataset of global solar irradiance (ghi), diffuse solar irradiance aka solar radiation bouncing of trees, clouds, etc (dhi), and cloud cover. I theorize that I can estimate the dhi given ghi and cloud cover for that day. Because more cloud cover could indicate more a portion of ghi from bouncing off clouds.

My goal is to build a prediction engine to predict DHI from GHI and Cloud Cover.

GHI values range from 0-1000 in and DHI from 0-500. Both are in w/m^2

Cloud Cover values range from 0 to 1.

A sample of the data looks as follows

Features:

ghi = [251 308 747 811 410 936 489 548 657  33  19  29 176  14 503  61 180   2
 487  62 283 156 224 704 771 187  99 696  17 810 426 113 205  22  14 710
  11 375  32 196 109   9 644 130 686 230 369 153 301  14  63  24  70 169
 191 525 576 177 384 271  60 648 551 594 186  45 850 117 507 133 791 192
 716 637  18  35 193 207 530 326  29 407  66  65 297 182 410 727 623  33
 178  12 820  32 424 925  39 681 473 451]

cloud_cover = [0.001 0.56  0.64  0.16  0.001 0.999 0.001 0.001 0.001 0.001 0.67  0.75
 0.75  0.66  0.001 0.84  0.08  0.999 0.44  0.001 0.001 0.001 0.09  0.35
 0.14  0.96  0.001 0.61  0.56  0.001 0.001 0.999 0.001 0.57  0.001 0.001
 0.999 0.001 0.41  0.44  0.15  0.96  0.999 0.001 0.999 0.001 0.9   0.999
 0.999 0.98  0.06  0.001 0.999 0.999 0.95  0.32  0.001 0.999 0.09  0.38
 0.87  0.08  0.001 0.001 0.01  0.001 0.999 0.75  0.28  0.999 0.001 0.73
 0.88  0.001 0.001 0.55  0.07  0.001 0.31  0.999 0.75  0.08  0.41  0.001
 0.97  0.999 0.001 0.001 0.11  0.001 0.001 0.25  0.61  0.81  0.7   0.999
 0.001 0.001 0.001 0.999]

response:
dhi = [ 64 211 132  96  79 142  88 155 133  21  19  24  54  14  84  61  44   2
 124  46 258  99 105 227 284 164  92 346  17 225  78 113  87  22  14 203
  10  55  25 188  40   9 460  84 433 135 152 111  95  14  34  21  43 167
 159 281 113 171 158 231  60 341  77 158  48  28  95  75  63  84 127 158
  74 137  14  25  92  56  68 294  21  63  40  49 268 163  72 174  78  27
 171  12 128  32  79 160  39 121 199 373]

And visually

I have fit a GLS regression with the following results

                                 GLS Regression Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.689
Model:                            GLS   Adj. R-squared (uncentered):              0.689
Method:                 Least Squares   F-statistic:                              4981.
Date:                Mon, 30 Mar 2020   Prob (F-statistic):                        0.00
Time:                        11:35:32   Log-Likelihood:                         -25872.
No. Observations:                4494   AIC:                                  5.175e+04
Df Residuals:                    4492   BIC:                                  5.176e+04
Df Model:                           2                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1             0.2163      0.003     73.546      0.000       0.211       0.222
x2            64.0145      2.227     28.741      0.000      59.648      68.381
==============================================================================
Omnibus:                      612.530   Durbin-Watson:                   0.526
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              994.138
Skew:                           0.933   Prob(JB):                    1.34e-216
Kurtosis:                       4.352   Cond. No.                         841.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

My major questions are:

1. How do I evaluate this model? The r2 is quite high, but when I run naive predictions, my results are on average way off from the actuals. Usually by 40%-50%.

2. How exactly do I specify the correct covariance matrix of errors? I am beginning to think this is part of the problem.

3. Is GLS even the best approach overall to counteract heteroscedasticity? I have tried nearly everything, including log and box-cox/power transformations. While they yield high r2 scores, the average error in predictions remains off by at least 30%.

Let me know if there is any more information I can provide. Thank you everyone who has helped me on this data journey!

EDIT:

The data is hourly measurements at one location for a period of one year. I have run regressions and p-tests for other variables including air pressure, but they do not appear to be statistically significant or they end up having lower r2 scores than without.

Answer 1

How where the measurements done? From where (ground or satelite, or combination)? Any other variables? (maybe time of day, ...) All measurements the same location, or different locations? We need much more detail, without that you will only be told generalities ...

Start out with plotting the data! In this case a conditioning plot may be informative:

(made in R with function coplot) The order of the panels are left to right, bottom up. Maybe the variance is increasing also with increasing cloud cover, or is this an artifact made by few data? Can you share a link to the full data? I really think you need to give more details.

And, to add, you should not plunge into modeling, start with a graphical analysis, using the full data. There might be much of interest to find! So, in a way I guess your question is premature.