Granger Causality Analysis of Solar Activity and Global Temperature: Linear and Nonlinear Tests on the Instrumental Record
Granger Causality Analysis of Solar Activity and Global Temperature: Linear and Nonlinear Tests on the Instrumental Record
Author: stepstep_labs
Abstract
We apply linear and nonlinear Granger causality testing to evaluate whether monthly sunspot numbers carry predictive information for global land-ocean temperature anomalies. Using the overlapping period of the SILSO v2.0 sunspot record and the NOAA NCEI temperature series, we analyze N = 2,100 monthly observations. Augmented Dickey-Fuller tests confirm that the sunspot series is stationary (ADF = −8.20, p < 0.01) while the temperature series contains a unit root (ADF = −1.69, p > 0.10); first-differencing renders the temperature series stationary (ADF = −19.86, p < 0.01). Information-criterion lag selection identifies p = 8 (AIC) and p = 6 (BIC) as optimal for the temperature autoregression. In level-form Granger tests across nine lag orders from 1 to 132 months, no specification achieves significance at the 5% level; the smallest p-value is 0.169 at lag 24. Granger tests on first-differenced series likewise yield no rejection, with the smallest p-value of 0.127 at the BIC-optimal lag of 5. A bidirectional falsification test confirms that temperature does not Granger-cause sunspot number (one isolated marginal rejection at p = 0.048 does not survive Bonferroni correction). Partial correlations between AR-filtered residuals are negligible (|r| < 0.013), and a 500-iteration permutation test at lag 12 returns a p-value of 0.513. Sliding-window analysis across twenty-six 50-year windows reveals only a single marginally significant window (1905–1954, p = 0.036), consistent with the expected false-positive rate. Three nonlinear extensions—squared-SSN Granger tests, epoch-split analyses, and regime-separated tests by solar-activity level—produce no robust evidence of nonlinear predictive content. These results indicate that solar cycle variability, as indexed by the international sunspot number, does not provide statistically detectable Granger-causal information for monthly global temperature anomalies over the instrumental record.
1. Introduction
The relationship between solar activity and Earth's climate has been debated for over a century. William Herschel proposed in 1801 that sunspot number might correlate with grain prices—a proxy for climate—and successive generations of researchers have sought to quantify the solar influence on temperature, precipitation, and atmospheric circulation. Total solar irradiance (TSI) varies by roughly 0.1% over the 11-year Schwabe cycle, corresponding to a radiative forcing of approximately 0.1–0.2 W m⁻² between cycle minimum and maximum. While detectable, this forcing is small relative to the approximately 2.7 W m⁻² attributed to anthropogenic greenhouse gas increases since the preindustrial era.
Spectral approaches have been the dominant methodology for investigating solar-climate linkages. Studies using Fourier, wavelet, and multitaper techniques have reported weak but sometimes statistically significant power near 11-year and 22-year (Hale cycle) periods in temperature proxies and instrumental records. However, spectral coherence at a particular frequency does not establish causal direction: both series could be driven by an unobserved common factor, or the apparent coherence could be an artifact of shared low-frequency trends.
Granger causality provides a complementary framework. Introduced by Granger (1969) and formalized by Sims (1972), the Granger test asks a sharply defined question: does the past of series X improve prediction of series Y beyond what the past of Y alone provides? This is a test of incremental predictive information, not of physical mechanism. Nonetheless, a failure to reject the null hypothesis of no Granger causality places a useful upper bound on the statistical detectability of any causal pathway, direct or indirect.
Several prior studies have applied Granger-type tests to solar-climate data, with mixed results. Triacca (2005) found no Granger causality from TSI to global temperature over the 20th century. Kaufmann and Stern (1997) reported weak evidence of solar influence that vanished when anthropogenic forcings were included. More recently, studies using nonlinear extensions of Granger causality (e.g., transfer entropy) have also failed to find robust solar-to-temperature information flow.
A methodological requirement for valid Granger testing is that the series be stationary or appropriately differenced to achieve stationarity. Global temperature anomalies are widely recognized as containing a unit root driven by the secular warming trend. Failure to account for non-stationarity can produce spurious Granger rejections, rendering standard F-tests unreliable. Additionally, the choice of lag order in Granger regressions must be justified through formal information criteria rather than ad hoc selection to avoid overfitting or underfitting the autoregressive structure.
This paper contributes to the literature by applying a comprehensive battery of linear and nonlinear Granger tests to the longest available monthly instrumental records: the SILSO v2.0 international sunspot number (1749–present) and the NOAA NCEI global land-ocean temperature anomaly series (1850–present). We address stationarity through Augmented Dickey-Fuller (ADF) testing and by conducting parallel analyses on both level-form and first-differenced series. Lag orders are selected via the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), with a broad multi-lag sensitivity analysis to ensure robustness. We extend beyond the linear framework through squared-SSN Granger tests, epoch-split analyses, and regime-separated tests stratified by solar-activity level. Additional robustness checks include bidirectional testing, permutation inference, partial correlations on AR-filtered residuals, sliding-window temporal stability analysis, and linear detrending.
2. Methods
2.1 Data
Sunspot number. We use the monthly mean total sunspot number from the World Data Center SILSO, version 2.0, maintained by the Royal Observatory of Belgium. The series is semicolon-delimited with columns for year, month, decimal year, monthly mean sunspot number (SSN), standard deviation, number of observations, and a definitive/provisional indicator. Over the common period, SSN has a mean of 83.74, a standard deviation of 67.79, and ranges from 0.0 to 359.4.
Temperature anomaly. We use the NOAA NCEI global (land and ocean combined) monthly temperature anomaly series, expressed in degrees Celsius relative to the 1901–2000 base period. Over the common period, the mean anomaly is +0.057 °C, the standard deviation is 0.368 °C, and values range from −0.67 °C to +1.40 °C.
Data were accessed from the SILSO and NOAA NCEI archives. The sunspot record spans 1749–present; the temperature record spans 1850–present. We use the overlapping period, yielding N = 2,100 paired monthly observations.
2.2 Unit Root Testing
Stationarity is a prerequisite for valid Granger causality inference. We apply the Augmented Dickey-Fuller (ADF) test to each series. The ADF regression is:
ΔY_t = α + γ Y_{t−1} + Σᵢ₌₁ᵏ δᵢ ΔY_{t−i} + εₜ
where k is the number of augmenting lags selected by AIC (searching over k = 0, 1, …, 24). The null hypothesis is H₀: γ = 0 (unit root, non-stationary). The test statistic is the t-ratio of γ̂, compared against Dickey-Fuller critical values (asymptotic values with constant: 1% = −3.43, 5% = −2.86, 10% = −2.57).
When non-stationarity is detected, we first-difference the affected series and confirm that the differenced series is stationary before proceeding with Granger tests on differenced data.
2.3 Lag Order Selection
Rather than selecting lag orders ad hoc, we determine the optimal autoregressive order for the restricted (AR-only) temperature model using the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC):
AIC(p) = n ln(RSS_p / n) + 2(p + 1)
BIC(p) = n ln(RSS_p / n) + (p + 1) ln(n)
where RSS_p is the residual sum of squares from the AR(p) model and n is the number of effective observations. We evaluate p = 1, 2, …, 36 and report the AIC- and BIC-optimal orders. Granger tests are conducted at the optimal lag orders, with additional lags spanning 1 to 132 months (approximately one full solar cycle) reported as a sensitivity analysis.
2.4 Granger Causality Test
For each lag order p, we estimate two models by ordinary least squares (OLS):
Restricted model (AR-only):
T_t = α₀ + Σᵢ₌₁ᵖ αᵢ T_{t−i} + εₜ
Unrestricted model (AR + lagged SSN):
T_t = α₀ + Σᵢ₌₁ᵖ αᵢ T_{t−i} + Σⱼ₌₁ᵖ βⱼ SSN_{t−j} + εₜ
The F-statistic for the joint null hypothesis H₀: β₁ = β₂ = … = βₚ = 0 is:
F = [(RSS_r − RSS_u) / p] / [RSS_u / (n − 2p − 1)]
where RSS_r and RSS_u denote the residual sum of squares from the restricted and unrestricted models. Under H₀, this statistic follows an F(p, n − 2p − 1) distribution. We conduct Granger tests in both level form and on first-differenced series.
OLS coefficients are computed via the normal equations: β̂ = (X'X)⁻¹X'y, with the matrix inverse obtained by Gauss-Jordan elimination with partial pivoting. P-values are computed from the regularized incomplete beta function representation of the F cumulative distribution function.
2.5 Bidirectional Testing
As a falsification check, we reverse the direction and test whether temperature anomaly Granger-causes sunspot number. Since there is no plausible physical mechanism by which Earth's surface temperature could influence solar magnetic activity, any significant result would indicate statistical artifact and call into question results in the forward direction.
2.6 Partial Correlation
For each AR order p ∈ {12, 24, 60, 132}, we fit separate AR(p) models to the sunspot and temperature series, extract the residuals, and compute the Pearson correlation between the two residual series. This measures the linear association between the "innovation" components of each series after their respective autoregressive structures have been removed.
2.7 Permutation Test
To validate the parametric F-test, we perform a permutation test at lag order p = 12. We compute the observed F-statistic, then randomly shuffle the sunspot series 500 times (preserving the temperature series intact), recomputing F for each permutation. The permutation p-value is (K + 1) / (B + 1), where K is the number of permutations yielding F ≥ F_observed and B = 500 is the number of permutations.
2.8 Sliding-Window Analysis
We compute the Granger test in rolling 50-year (600-month) windows, stepping forward by 5 years (60 months), yielding 26 windows. This is done at two lag orders: p = 12 (one year) and p = 60 (five years). The purpose is to detect whether the solar-temperature relationship is temporally unstable, potentially masked by epoch-specific confounders.
2.9 Detrended Analysis
To address the concern that the monotonic warming trend in the temperature series could obscure or create spurious correlations, we remove a least-squares linear trend from the temperature series and repeat the Granger tests.
2.10 Nonlinear Extensions
Linear Granger causality may fail to capture complex, nonlinear interactions between solar cycles and climate systems. We implement three nonlinear extensions:
Squared-SSN Granger test. We replace the SSN series with SSN² in the unrestricted Granger regression, testing whether the magnitude (rather than the level) of solar activity carries predictive content for temperature.
Epoch-split analysis. We divide the record at the midpoint into two equal halves, independently detrend each half, and run Granger tests separately. This allows the solar-temperature relationship to differ across historical periods with distinct forcing regimes.
Regime-separated Granger test. Using a 132-month (one solar cycle) rolling average of SSN, we classify months into "high-activity" and "low-activity" regimes relative to the median. Granger tests are then conducted separately within contiguous regime blocks exceeding 200 months, after independent detrending within each block. This captures potential threshold effects whereby solar influence on climate operates differently during active versus quiet solar epochs.
Additionally, we compute the Pearson correlation between the absolute values of AR-filtered residuals from both series, testing whether the volatility (rather than level) of solar innovations is associated with the volatility of temperature innovations.
2.11 Cross-Correlation Function
We compute the sample cross-correlation function (CCF) between SSN and temperature anomaly at integer lags from −240 to +240 months (±20 years). Positive lag k means sunspot number leads temperature by k months.
3. Results
3.1 Unit Root Tests
Table 1 reports the ADF test results for each series.
Table 1. Augmented Dickey-Fuller unit root tests
| Series | ADF statistic | Optimal lag (AIC) | 1% critical | 5% critical | 10% critical | Conclusion |
|---|---|---|---|---|---|---|
| SSN (level) | −8.197 | 24 | −3.43 | −2.86 | −2.57 | Stationary (p < 0.01) |
| Temperature (level) | −1.693 | 7 | −3.43 | −2.86 | −2.57 | Non-stationary |
| ΔTemperature | −19.859 | 10 | −3.43 | −2.86 | −2.57 | Stationary (p < 0.01) |
| ΔSSN | −6.986 | 24 | −3.43 | −2.86 | −2.57 | Stationary (p < 0.01) |
The sunspot number is stationary in levels, consistent with its quasi-periodic oscillatory behavior. The temperature anomaly series fails to reject the unit root null at even the 10% level (ADF = −1.693), confirming the well-known non-stationarity driven by the secular warming trend. First-differencing the temperature series yields strong rejection of the unit root (ADF = −19.859), confirming that the series is integrated of order one, I(1). The first-differenced SSN is also stationary (ADF = −6.986).
These results motivate a dual analytical strategy: we conduct Granger tests on both level-form series (for comparability with earlier literature) and on first-differenced series (to satisfy the stationarity requirement for valid inference).
3.2 Information-Criterion Lag Selection
Table 2 reports AIC and BIC values for the temperature AR(p) model across selected lag orders.
Table 2. Information criteria for temperature AR(p) model
| Lag order (p) | AIC | BIC |
|---|---|---|
| 1 | −9,186.0 | −9,174.7 |
| 3 | −9,497.0 | −9,474.4 |
| 6 | −9,541.4 | −9,501.9 |
| 8 | −9,546.1 | −9,495.3 |
| 12 | −9,545.2 | −9,471.8 |
| 24 | −9,518.9 | −9,377.9 |
| 36 | −9,489.2 | −9,280.8 |
AIC selects p = 8 and BIC selects p = 6 as the optimal lag orders for the temperature autoregression. Both criteria penalize models beyond approximately p = 8, indicating that longer lag structures do not meaningfully improve fit after adjusting for model complexity. For the first-differenced temperature series, AIC selects p = 11 and BIC selects p = 5.
3.3 Granger Causality: Sunspot → Temperature (Levels)
Table 3 presents the F-statistics and p-values for the primary Granger test (SSN → Temperature) in levels, at information-criterion-optimal and sensitivity-analysis lag orders.
Table 3. Granger causality test: SSN → Temperature (levels)
| Lag order (p) | F-statistic | p-value | RSS (restricted) | RSS (unrestricted) | n | Note |
|---|---|---|---|---|---|---|
| 1 | 0.017 | 0.898 | 26.336 | 26.336 | 2,099 | |
| 6 | 1.461 | 0.188 | 21.837 | 21.745 | 2,094 | BIC-optimal |
| 8 | 1.447 | 0.172 | 21.631 | 21.511 | 2,092 | AIC-optimal |
| 12 | 0.943 | 0.502 | 21.329 | 21.213 | 2,088 | |
| 24 | 1.274 | 0.169 | 20.675 | 20.368 | 2,076 | |
| 36 | 0.913 | 0.618 | 20.067 | 19.741 | 2,064 | |
| 60 | 0.858 | 0.773 | 19.701 | 19.187 | 2,040 | |
| 120 | 0.984 | 0.533 | 18.230 | 17.071 | 1,980 | |
| 132 | 1.006 | 0.467 | 17.931 | 16.634 | 1,968 |
No lag order achieves significance at the 5% level. At the BIC-optimal lag of 6, F = 1.461 and p = 0.188. At the AIC-optimal lag of 8, F = 1.447 and p = 0.172. The smallest p-value across all lags is 0.169 at p = 24. At the solar-cycle lag of p = 132, F = 1.006 and p = 0.467, indicating no detectable predictive content. The F-statistics cluster around 1.0 across all lag orders, precisely what is expected under the null hypothesis.
3.4 Granger Causality: First-Differenced Series
Table 4 presents Granger test results on first-differenced series (ΔSSN → ΔTemperature), which satisfy the stationarity requirement.
Table 4. Granger causality test: ΔSSN → ΔTemperature (first-differenced)
| Lag order (p) | F-statistic | p-value | n | Note |
|---|---|---|---|---|
| 1 | 0.043 | 0.835 | 2,098 | |
| 5 | 1.719 | 0.127 | 2,094 | BIC-optimal for ΔTemp |
| 6 | 1.588 | 0.147 | 2,093 | |
| 11 | 1.010 | 0.435 | 2,088 | AIC-optimal for ΔTemp |
| 12 | 1.104 | 0.352 | 2,087 | |
| 24 | 1.277 | 0.166 | 2,075 | |
| 36 | 0.949 | 0.557 | 2,063 | |
| 60 | 0.927 | 0.636 | 2,039 |
No lag order achieves significance at the 5% level. The smallest p-value is 0.127 at the BIC-optimal lag of 5. This confirms that the null result from the level-form analysis is not an artifact of non-stationarity: even after differencing both series to achieve stationarity, sunspot changes carry no detectable predictive information for temperature changes.
3.5 Bidirectional Test: Temperature → Sunspot
Table 5 presents the reverse-direction Granger test.
Table 5. Granger causality test: Temperature → SSN (levels)
| Lag order (p) | F-statistic | p-value |
|---|---|---|
| 1 | 0.030 | 0.862 |
| 6 | 2.119 | 0.048 |
| 8 | 1.541 | 0.138 |
| 12 | 1.053 | 0.397 |
| 24 | 0.721 | 0.835 |
| 36 | 0.883 | 0.669 |
| 60 | 1.031 | 0.411 |
| 120 | 0.961 | 0.603 |
| 132 | 0.961 | 0.609 |
Only one lag order (p = 6) achieves nominal significance at the 5% level, with F = 2.119 and p = 0.048. Given that nine tests are conducted, the Bonferroni-corrected threshold is 0.05/9 = 0.0056, and this result does not survive correction. The isolated rejection is consistent with the expected false-positive rate. Importantly, temperature cannot physically influence solar magnetic activity, so this single marginal result confirms the tests are well-behaved and not systematically anti-conservative.
3.6 Partial Correlation
After removing autoregressive structure from both series, the residual correlations are negligible:
| AR order | Partial r | t-statistic | n |
|---|---|---|---|
| 12 | 0.009 | 0.428 | 2,088 |
| 24 | 0.012 | 0.553 | 2,076 |
| 60 | 0.005 | 0.214 | 2,040 |
| 132 | 0.002 | 0.104 | 1,968 |
All partial correlations are below 0.013 in absolute value, and no t-statistic exceeds 0.55. The innovation components of the two series are uncorrelated.
3.7 Permutation Test
At lag order p = 12, the observed F-statistic is 0.943. Among 500 random permutations of the sunspot series, 256 yielded an F-statistic equal to or exceeding the observed value, giving a permutation p-value of 0.513. This is consistent with the parametric p-value of 0.502 and confirms that the null hypothesis cannot be rejected.
3.8 Sliding-Window Analysis
At lag order p = 12, the 26 windows show F-statistics ranging from 0.600 (1880–1929) to 1.863 (1905–1954). Only the 1905–1954 window achieves nominal significance (p = 0.036), but this single rejection among 26 tests is expected by chance alone (expected false positives at α = 0.05: 26 × 0.05 = 1.3). No consistent temporal pattern of significance is evident—windows before, during, and after the mid-20th century all show non-significant results.
At lag order p = 60, no window achieves significance. The F-statistics range from 0.731 (1965–2014) to 1.079 (1945–1994), with all p-values exceeding 0.33. The relationship is uniformly null across the entire span.
3.9 Detrended Analysis
After removing the linear warming trend (0.056 °C per decade), the contemporaneous Pearson correlation between SSN and detrended temperature is r = −0.062. The Granger test results on the detrended series are substantively identical to the level-form results:
Table 6. Granger causality test on detrended temperature
| Lag order (p) | F-statistic | p-value |
|---|---|---|
| 1 | 0.301 | 0.584 |
| 6 | 1.408 | 0.208 |
| 8 | 1.403 | 0.190 |
| 12 | 0.925 | 0.521 |
| 24 | 1.266 | 0.175 |
| 36 | 0.909 | 0.625 |
| 60 | 0.853 | 0.783 |
| 120 | 0.974 | 0.564 |
| 132 | 0.997 | 0.495 |
No lag order is significant. Detrending slightly reduces the F-statistics at most lag orders, confirming that the warming trend was not creating spurious Granger causality.
3.10 Cross-Correlation Function
The contemporaneous cross-correlation between SSN and temperature anomaly is effectively zero (r = −0.003 at lag 0). Over the full range of lags examined (±240 months), the peak positive value is r = 0.171 at lag +240 months (sunspot leading temperature by 20 years) and the peak negative value is r = −0.033 at lag +34 months. At shorter lags (−60 to +60 months), the maximum absolute CCF is 0.033. The increasing CCF at large positive lags is an artifact of the secular warming trend in the temperature record, not a genuine lagged solar influence.
3.11 Nonlinear Extensions
Squared-SSN Granger test. Table 7 presents results for testing whether SSN² Granger-causes temperature.
Table 7. Granger causality test: SSN² → Temperature
| Lag order (p) | F-statistic | p-value |
|---|---|---|
| 1 | 0.280 | 0.597 |
| 6 | 2.347 | 0.029 |
| 12 | 1.312 | 0.204 |
| 24 | 1.388 | 0.100 |
At p = 6, the squared-SSN test yields a nominally significant result (F = 2.347, p = 0.029). However, this single result among four tests does not survive Bonferroni correction (threshold = 0.0125), and significance vanishes at all other lag orders. The isolated rejection at p = 6 is consistent with the expected false-positive rate across the full suite of nonlinear tests.
Epoch-split analysis. Splitting the record at the midpoint and independently detrending each half, no Granger test achieves significance in either the early or the late epoch:
| Epoch | Lag 6 (p-value) | Lag 12 (p-value) | Lag 24 (p-value) |
|---|---|---|---|
| First half (1850–1937) | 0.510 | 0.746 | 0.475 |
| Second half (1937–present) | 0.321 | 0.277 | 0.373 |
The null result is consistent across distinct historical periods with different volcanic, ENSO, and anthropogenic forcing regimes.
Regime-separated Granger test. Classifying months into high- and low-activity regimes based on the 132-month rolling SSN (median = 88.5), we identify four contiguous blocks exceeding 200 months: three high-activity epochs (1860–1880, 1941–1971, 1979–2009) and one low-activity epoch (1880–1940). Table 8 reports results at lag 12 within each regime, after independent detrending.
Table 8. Regime-separated Granger causality test (SSN → Temperature, lag 12, detrended)
| Regime | Period | n | F-statistic | p-value |
|---|---|---|---|---|
| High activity | 1860–1880 | 232 | 0.777 | 0.674 |
| High activity | 1941–1971 | 360 | 0.798 | 0.652 |
| High activity | 1979–2009 | 357 | 1.227 | 0.263 |
| Low activity | 1880–1940 | 727 | 0.666 | 0.785 |
No regime yields a significant result. The null finding is invariant to whether solar activity is high or low, providing no evidence for threshold or asymmetric nonlinear effects.
Absolute-residual correlation. After AR(12) filtering, the correlation between |SSN residuals| and |temperature residuals| is r = −0.013 (t = −0.601, n = 2,088); after AR(24) filtering, it is r = −0.005 (t = −0.218, n = 2,076). The volatility components of the two series are uncorrelated, ruling out variance-level nonlinear coupling.
4. Discussion
4.1 Interpretation of the Null Result
The central finding of this study is unambiguous: the international sunspot number does not Granger-cause global monthly temperature anomalies over the instrumental record. This null result is robust to the choice of lag order (information-criterion-optimal or sensitivity range from 1 to 132 months), to the treatment of non-stationarity (level-form, first-differenced, and detrended series), to the functional form of the predictor (linear SSN, squared SSN, regime-separated), to permutation-based inference, and to temporal subsampling via sliding windows.
This finding does not imply that the Sun has zero influence on climate. It implies that the signal, if present, is too small to be detected by Granger-type frameworks given the noise properties of the monthly temperature series. The temperature series is dominated by strong serial autocorrelation, internal climate variability (ENSO, volcanic aerosols, Atlantic Multidecadal Oscillation), and the monotonic anthropogenic warming trend. Against this background, the ~0.1% solar irradiance variation translates to a radiative forcing perturbation that is an order of magnitude smaller than the dominant anthropogenic signal.
4.2 The Importance of Stationarity
The ADF tests reveal a critical asymmetry: the sunspot series is stationary, while the temperature series contains a unit root. This asymmetry is consequential for Granger inference. If both series were I(1) and cointegrated, error-correction models would be appropriate; if both were stationary, standard F-tests would suffice. In the mixed-integration case, standard Granger tests applied to levels can be unreliable—they may reject the null spuriously due to the stochastic trend in temperature, or they may fail to reject because the nonstationary component dominates the regression.
Our dual strategy resolves this ambiguity. The level-form tests (Table 3) yield no rejections, but their validity is qualified by the unit root in temperature. The first-differenced tests (Table 4) are conducted on stationary series and are therefore fully valid under standard asymptotic theory. The fact that both approaches yield the same null result—with remarkably similar p-values—provides strong evidence that the finding is genuine and not an artifact of non-stationarity.
4.3 Lag Selection and Sensitivity
Information criteria identify p = 6 (BIC) and p = 8 (AIC) as the optimal lag orders for the temperature autoregression, corresponding to 6–8 months of memory in the temperature process. These short optimal lags reflect the dominance of month-to-month persistence in the temperature record, with diminishing returns beyond approximately half a year. At these optimal lags, the Granger p-values are 0.188 (BIC) and 0.172 (AIC), both far from significance.
The multi-lag sensitivity analysis from p = 1 to p = 132 serves two purposes: it checks whether significance might emerge at lags beyond the information-criterion optimum, and it allows for the possibility that the solar signal operates on longer timescales (e.g., the full 11-year cycle) that AIC/BIC penalize due to the large number of additional parameters. The absence of significance at any lag order, including p = 132, rules out this possibility.
4.4 Comparison with Spectral Studies
The spectral literature has reported detectable solar-cycle power in some climate records, particularly in pre-industrial proxies and stratospheric temperatures. Our results are not contradictory. Spectral methods detect periodic power, which may be present even when predictive information (in the Granger sense) is negligible. If the solar signal contributes 0.05–0.1 °C of periodic variation against monthly noise with a standard deviation of 0.37 °C, the signal-to-noise ratio is below 0.3, making Granger detection infeasible at monthly resolution.
4.5 Nonlinear Extensions
The three nonlinear tests collectively address the concern that linear Granger causality may miss complex solar-climate interactions. The squared-SSN test captures intensity-dependent effects; the epoch-split analysis allows for time-varying coupling coefficients; and the regime-separated test permits threshold behavior where solar influence operates differently during active versus quiet epochs. None of these extensions produces robust evidence of predictive content. The single nominally significant result (SSN² at lag 6, p = 0.029) does not survive multiple-testing correction and is inconsistent with the null results at all other lag orders and functional forms. The absolute-residual correlation analysis further rules out variance-level coupling. While these nonlinear extensions do not exhaust the space of possible nonlinear specifications (e.g., neural network Granger tests, transfer entropy), they address the most physically motivated nonlinear hypotheses regarding solar-climate interactions.
4.6 Sliding-Window Stability
The sliding-window analysis reveals that the null result is temporally stable. The single marginal rejection in the 1905–1954 window (p = 0.036 at lag 12) does not constitute credible evidence for an epoch-specific solar influence. With 26 independent tests at α = 0.05, the expected number of false positives is 1.3, and one rejection is squarely within this expectation. No clustering of near-significant results is observed in any sub-period.
4.7 Limitations
Several limitations should be acknowledged. First, although we implemented three nonlinear extensions, the space of possible nonlinear specifications is vast. Transfer entropy or neural-network-based Granger tests could be explored in future work. Second, we use the sunspot number as a proxy for solar forcing; direct TSI reconstructions, available from 1978 onward, might yield different results over the satellite era, though the shorter record would substantially reduce statistical power. Third, the Granger framework tests pairwise causality and cannot account for confounders that affect both series simultaneously—volcanic aerosol loading, for example, could mask solar signals if eruptions coincide with particular solar cycle phases. Fourth, the analysis is conducted at monthly resolution; aggregation to annual or decadal scales might reveal signals obscured by high-frequency climate noise, though at the cost of drastically reduced sample sizes. Fifth, the international sunspot number has undergone several recalibrations, most recently the SILSO v2.0 update, and any remaining inhomogeneities could affect results, though the recalibration primarily shifted the absolute scale rather than the temporal pattern.
5. Conclusion
We find no evidence that the international sunspot number Granger-causes global monthly temperature anomalies. Augmented Dickey-Fuller tests confirm that the temperature series is non-stationary (ADF = −1.69) while the sunspot series is stationary (ADF = −8.20); first-differencing renders the temperature series stationary (ADF = −19.86). Information criteria select lag orders of 6 (BIC) and 8 (AIC) for the temperature autoregression. The null hypothesis of no predictive information from sunspots to temperature cannot be rejected at any lag order from 1 to 132 months, whether applied to level-form, first-differenced, or detrended temperature series, and whether tested across the full record or within 50-year subperiods. Three nonlinear extensions—squared-SSN Granger tests, epoch-split analyses, and regime-separated tests—produce no robust evidence of nonlinear predictive content. Partial correlations between AR-filtered innovations are negligible (|r| < 0.013), and a permutation test independently confirms the null. The bidirectional falsification test behaves as expected.
These results are consistent with the scientific consensus that anthropogenic greenhouse gas forcing has overwhelmingly dominated the global temperature trend over the instrumental period, rendering the ~0.1% solar irradiance cycle undetectable by time-series predictive tests at monthly resolution. Solar variability may still contribute a small periodic component to climate, but this component is not large enough to improve upon a purely autoregressive temperature forecast, even when the model is given 11 years of lagged sunspot information or when nonlinear functional forms are considered.
References
Granger, C. W. J. (1969). "Investigating Causal Relations by Econometric Models and Cross-spectral Methods." Econometrica, 37(3), 424–438.
Sims, C. A. (1972). "Money, Income, and Causality." American Economic Review, 62(4), 540–552.
Triacca, U. (2005). "Is Granger causality analysis appropriate to investigate the relationship between atmospheric concentration of carbon dioxide and global surface air temperature?" Theoretical and Applied Climatology, 81, 133–135.
Kaufmann, R. K., & Stern, D. I. (1997). "Evidence for human influence on climate from hemispheric temperature relations." Nature, 388, 39–44.
Dickey, D. A., & Fuller, W. A. (1979). "Distribution of the Estimators for Autoregressive Time Series with a Unit Root." Journal of the American Statistical Association, 74(366), 427–431.
Akaike, H. (1974). "A New Look at the Statistical Model Identification." IEEE Transactions on Automatic Control, 19(6), 716–723.
Schwarz, G. (1978). "Estimating the Dimension of a Model." Annals of Statistics, 6(2), 461–464.
SILSO World Data Center. "International Sunspot Number Monthly Bulletin and Online Catalogue." Royal Observatory of Belgium. Retrieved from https://www.sidc.be/SILSO/
NOAA National Centers for Environmental Information. "Global Surface Temperature Anomalies." Retrieved from https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/
Lean, J., & Rind, D. (2008). "How natural and anthropogenic influences alter global and regional surface temperatures: 1889 to 2006." Geophysical Research Letters, 35, L18701.
Gray, L. J., et al. (2010). "Solar influences on climate." Reviews of Geophysics, 48, RG4001.
Herschel, W. (1801). "Observations Tending to Investigate the Nature of the Sun." Philosophical Transactions of the Royal Society of London, 91, 265–318.
Friis-Christensen, E., & Lassen, K. (1991). "Length of the Solar Cycle: An Indicator of Solar Activity Closely Associated with Climate." Science, 254(5032), 698–700.
Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
# Skill: Granger Causality — Sunspot Number → Global Temperature Anomaly (with ADF, AIC/BIC, Nonlinear Extensions)
## allowed-tools
Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(echo *)
## Overview
This skill tests whether the international sunspot number (SILSO v2.0) Granger-causes global land-ocean temperature anomalies (NOAA NCEI) over the instrumental record. The analysis includes Augmented Dickey-Fuller unit root tests, AIC/BIC-based lag selection, Granger tests on both level-form and first-differenced series, bidirectional falsification, partial correlation on AR residuals, a permutation test, sliding-window stability analysis, detrended robustness checks, and nonlinear extensions (squared-SSN, epoch-split, and regime-separated Granger tests).
## Data Requirements
1. **Sunspot data**: SILSO v2.0 monthly sunspot number CSV (semicolon-delimited; columns: year;month;decimal_year;ssn;...). Download from https://www.sidc.be/SILSO/INFO/snmtotcsv.php
2. **Temperature data**: NOAA NCEI global land-ocean monthly temperature anomaly. Download from https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land_ocean/1/0/1850-2024/data.csv
## Dependencies
- Pure Python standard library only (math, random, csv, os)
- random.seed(42) for reproducibility
## Analysis Script
```python
#!/usr/bin/env python3
"""
Granger Causality Analysis: Solar Activity → Global Temperature Anomaly
Includes ADF tests, AIC/BIC lag selection, differenced-series Granger,
and nonlinear extensions.
Pure Python stdlib implementation.
"""
import math
import random
import csv
import os
random.seed(42)
# ─────────────────────────────────────────────────────────
# 1. LINEAR ALGEBRA HELPERS
# ─────────────────────────────────────────────────────────
def mat_inv(M):
"""Gauss-Jordan inverse with partial pivoting."""
n = len(M)
aug = [M[i][:] + [1.0 if j == i else 0.0 for j in range(n)] for i in range(n)]
for col in range(n):
mx = abs(aug[col][col]); mi = col
for row in range(col+1, n):
if abs(aug[row][col]) > mx:
mx = abs(aug[row][col]); mi = row
aug[col], aug[mi] = aug[mi], aug[col]
piv = aug[col][col]
if abs(piv) < 1e-15:
raise ValueError("Singular matrix")
for j in range(2*n):
aug[col][j] /= piv
for row in range(n):
if row == col:
continue
f = aug[row][col]
for j in range(2*n):
aug[row][j] -= f * aug[col][j]
return [aug[i][n:] for i in range(n)]
def ols_fit(X, y):
"""OLS via normal equations. Returns (beta, RSS)."""
n = len(y)
k = len(X[0])
XtX = [[0.0]*k for _ in range(k)]
Xty = [0.0]*k
for i in range(n):
xi = X[i]
yi = y[i]
for a in range(k):
xia = xi[a]
Xty[a] += xia * yi
for b in range(a, k):
v = xia * xi[b]
XtX[a][b] += v
if a != b:
XtX[b][a] += v
XtX_inv = mat_inv(XtX)
beta = [0.0]*k
for i in range(k):
s = 0.0
for j in range(k):
s += XtX_inv[i][j] * Xty[j]
beta[i] = s
rss = 0.0
for i in range(n):
yhat = sum(xi_j * beta[j] for j, xi_j in enumerate(X[i]))
rss += (y[i] - yhat)**2
return beta, rss
# ─────────────────────────────────────────────────────────
# 2. STATISTICAL HELPERS
# ─────────────────────────────────────────────────────────
def _log_gamma(x):
if x <= 0:
return float('inf')
coefs = [76.18009172947146, -86.50532032941677, 24.01409824083091,
-1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5]
tmp = x + 5.5
ser = 1.000000000190015
for j, c in enumerate(coefs):
ser += c / (x + 1 + j)
return (x + 0.5) * math.log(tmp) - tmp + math.log(2.5066282746310005 * ser / x)
def _beta_cf(a, b, x, max_iter=200, eps=1e-12):
qab = a + b; qap = a + 1.0; qam = a - 1.0
c = 1.0
d = 1.0 - qab * x / qap
if abs(d) < 1e-30: d = 1e-30
d = 1.0 / d; h = d
for m in range(1, max_iter + 1):
m2 = 2 * m
aa = m * (b - m) * x / ((qam + m2) * (a + m2))
d = 1.0 + aa * d
if abs(d) < 1e-30: d = 1e-30
c = 1.0 + aa / c
if abs(c) < 1e-30: c = 1e-30
d = 1.0 / d; h *= d * c
aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
d = 1.0 + aa * d
if abs(d) < 1e-30: d = 1e-30
c = 1.0 + aa / c
if abs(c) < 1e-30: c = 1e-30
d = 1.0 / d; delta = d * c; h *= delta
if abs(delta - 1.0) < eps:
return h
return h
def _betainc(a, b, x):
if x < 0.0 or x > 1.0: return 0.0
if x == 0.0 or x == 1.0: return x
lbeta = _log_gamma(a + b) - _log_gamma(a) - _log_gamma(b)
front = math.exp(lbeta + a * math.log(x) + b * math.log(1.0 - x))
if x < (a + 1.0) / (a + b + 2.0):
return front * _beta_cf(a, b, x) / a
else:
return 1.0 - front * _beta_cf(b, a, 1.0 - x) / b
def f_pvalue(f_val, d1, d2):
if f_val <= 0: return 1.0
x = d1 * f_val / (d1 * f_val + d2)
return 1.0 - _betainc(d1 / 2.0, d2 / 2.0, x)
def mean(x):
return sum(x) / len(x)
def std(x):
m = mean(x)
return math.sqrt(sum((xi - m)**2 for xi in x) / (len(x) - 1))
def pearson_r(x, y):
n = len(x); mx, my = mean(x), mean(y); sx, sy = std(x), std(y)
if sx == 0 or sy == 0: return 0.0
return sum((x[i]-mx)*(y[i]-my) for i in range(n)) / ((n-1)*sx*sy)
# ─────────────────────────────────────────────────────────
# 3. DATA LOADING
# ─────────────────────────────────────────────────────────
BASE = os.path.dirname(os.path.abspath(__file__))
def load_sunspot(path):
data = {}
with open(path, 'r') as f:
for line in f:
line = line.strip()
if not line or line.startswith('#'): continue
parts = line.split(';')
if len(parts) < 4: continue
try:
yr = int(parts[0].strip()); mo = int(parts[1].strip())
ssn = float(parts[3].strip())
except (ValueError, IndexError): continue
if ssn < 0: ssn = 0.0
data[(yr, mo)] = ssn
return data
def load_temperature(path):
data = {}
with open(path, 'r') as f:
for line in f:
line = line.strip()
if not line or line.startswith('#') or line.startswith('Date'): continue
parts = line.split(',')
if len(parts) < 2: continue
try:
date_str = parts[0].strip()
anom = float(parts[1].strip())
yr = int(date_str[:4]); mo = int(date_str[4:])
except (ValueError, IndexError): continue
data[(yr, mo)] = anom
return data
def load_gistemp(path):
data = {}
with open(path, 'r') as f:
lines = f.readlines()
header_idx = None
for i, line in enumerate(lines):
if line.strip().startswith('Year'):
header_idx = i; break
if header_idx is None: return data
for line in lines[header_idx+1:]:
parts = line.strip().split(',')
if len(parts) < 13: continue
try: yr = int(parts[0])
except ValueError: continue
for mo in range(1, 13):
val = parts[mo].strip()
if val == '***' or val == '': continue
try: data[(yr, mo)] = float(val)
except ValueError: continue
return data
# Load data
sunspot_path = os.path.join(os.path.dirname(BASE), 'sunspot_cusum', 'monthly_sunspot.csv')
if not os.path.exists(sunspot_path):
sunspot_path = os.path.join(BASE, '..', 'sunspot_cusum', 'monthly_sunspot.csv')
ssn_data = load_sunspot(sunspot_path)
noaa_path = os.path.join(BASE, 'temperature_data.csv')
gistemp_path = os.path.join(BASE, 'gistemp.csv')
temp_data = load_temperature(noaa_path)
temp_source = "NOAA NCEI"
if len(temp_data) < 100:
temp_data = load_gistemp(gistemp_path)
temp_source = "NASA GISTEMP v4"
common_keys = sorted(set(ssn_data.keys()) & set(temp_data.keys()))
ssn_series = [ssn_data[k] for k in common_keys]
temp_series = [temp_data[k] for k in common_keys]
N = len(ssn_series)
print(f"Data: {N} observations, {common_keys[0][0]}–{common_keys[-1][0]}")
print(f"Sources: SILSO v2.0, {temp_source}")
# ─────────────────────────────────────────────────────────
# 4. AUGMENTED DICKEY-FULLER TEST
# ─────────────────────────────────────────────────────────
def adf_test(series, max_lags=24, name="series"):
"""ADF test with AIC-selected augmentation lag."""
n = len(series)
dy = [series[i] - series[i-1] for i in range(1, n)]
best_aic = float('inf'); best_lag = 0; best_result = None
for p_lag in range(0, max_lags + 1):
y_dep = []; X_mat = []
for t in range(p_lag, len(dy)):
y_dep.append(dy[t])
row = [1.0, series[t]] # constant + Y_{t-1}
for i in range(1, p_lag + 1):
row.append(dy[t - i])
X_mat.append(row)
n_obs = len(y_dep); k = len(X_mat[0])
if n_obs <= k + 1: continue
try: beta, rss = ols_fit(X_mat, y_dep)
except: continue
aic = n_obs * math.log(rss / n_obs) + 2 * k
if aic < best_aic:
best_aic = aic; best_lag = p_lag
sigma2 = rss / (n_obs - k)
XtX = [[0.0]*k for _ in range(k)]
for i in range(n_obs):
xi = X_mat[i]
for a in range(k):
for b in range(a, k):
v = xi[a] * xi[b]
XtX[a][b] += v
if a != b: XtX[b][a] += v
try:
XtX_inv = mat_inv(XtX)
se = math.sqrt(sigma2 * XtX_inv[1][1])
best_result = (beta[1] / se, best_lag, n_obs)
except: continue
crit = {1: -3.43, 5: -2.86, 10: -2.57}
if best_result:
adf_stat, lag, n_obs = best_result
reject = "STATIONARY" if adf_stat < crit[5] else "NON-STATIONARY"
print(f" ADF({name}): stat={adf_stat:.4f}, lag={lag}, → {reject}")
return (adf_stat, lag, n_obs, reject)
return None
print("\n--- Unit Root Tests ---")
adf_ssn = adf_test(ssn_series, name="SSN")
adf_temp = adf_test(temp_series, name="Temp")
temp_diff = [temp_series[i] - temp_series[i-1] for i in range(1, N)]
ssn_diff = [ssn_series[i] - ssn_series[i-1] for i in range(1, N)]
adf_dtemp = adf_test(temp_diff, name="ΔTemp")
adf_dssn = adf_test(ssn_diff, name="ΔSSN")
# ─────────────────────────────────────────────────────────
# 5. AIC/BIC LAG SELECTION
# ─────────────────────────────────────────────────────────
def select_lag(y_series, max_lag=36):
"""Select AR lag by AIC and BIC."""
n = len(y_series); results = []
for p in range(1, max_lag + 1):
y_dep = []; X_mat = []
for t in range(p, n):
y_dep.append(y_series[t])
row = [1.0] + [y_series[t-i] for i in range(1, p+1)]
X_mat.append(row)
n_obs = len(y_dep); k = p + 1
if n_obs <= k + 1: continue
try: _, rss = ols_fit(X_mat, y_dep)
except: continue
aic = n_obs * math.log(rss / n_obs) + 2 * k
bic = n_obs * math.log(rss / n_obs) + k * math.log(n_obs)
results.append((p, aic, bic))
best_aic = min(results, key=lambda x: x[1])
best_bic = min(results, key=lambda x: x[2])
return best_aic[0], best_bic[0]
aic_p, bic_p = select_lag(temp_series)
print(f"\nLag selection (Temp): AIC={aic_p}, BIC={bic_p}")
aic_p_d, bic_p_d = select_lag(temp_diff)
print(f"Lag selection (ΔTemp): AIC={aic_p_d}, BIC={bic_p_d}")
# ─────────────────────────────────────────────────────────
# 6. GRANGER CAUSALITY TEST
# ─────────────────────────────────────────────────────────
def granger_test(y_series, x_series, p):
n = len(y_series)
if p >= n // 2: return None
y_dep = []; X_r = []; X_u = []
for t in range(p, n):
y_dep.append(y_series[t])
row_r = [1.0] + [y_series[t-i] for i in range(1, p+1)]
X_r.append(row_r)
row_u = row_r[:] + [x_series[t-j] for j in range(1, p+1)]
X_u.append(row_u)
n_obs = len(y_dep)
try:
_, rss_r = ols_fit(X_r, y_dep)
_, rss_u = ols_fit(X_u, y_dep)
except: return None
dof_den = n_obs - 2*p - 1
if dof_den <= 0 or rss_u <= 0: return None
F = ((rss_r - rss_u) / p) / (rss_u / dof_den)
return (F, f_pvalue(F, p, dof_den), rss_r, rss_u, n_obs)
# Primary Granger tests (levels)
all_lags = sorted(set([aic_p, bic_p, 1, 6, 12, 24, 36, 60, 120, 132]))
print(f"\n--- Granger: SSN → Temp (levels) ---")
for p in all_lags:
res = granger_test(temp_series, ssn_series, p)
if res:
F, pval, rss_r, rss_u, n_obs = res
note = ""
if p == aic_p: note += " [AIC]"
if p == bic_p: note += " [BIC]"
sig = "*" if pval < 0.05 else ""
print(f" p={p:>3d}: F={F:.4f}, p={pval:.6f} {sig}{note}")
# Granger tests on differenced series
diff_lags = sorted(set([aic_p_d, bic_p_d, 1, 6, 12, 24, 36, 60]))
print(f"\n--- Granger: ΔSSN → ΔTemp (differenced) ---")
for p in diff_lags:
res = granger_test(temp_diff, ssn_diff, p)
if res:
F, pval, _, _, n_obs = res
note = ""
if p == aic_p_d: note += " [AIC]"
if p == bic_p_d: note += " [BIC]"
sig = "*" if pval < 0.05 else ""
print(f" p={p:>3d}: F={F:.4f}, p={pval:.6f} {sig}{note}")
# Bidirectional test
print(f"\n--- Bidirectional: Temp → SSN (levels) ---")
for p in all_lags:
res = granger_test(ssn_series, temp_series, p)
if res:
F, pval, _, _, _ = res
sig = "*" if pval < 0.05 else ""
print(f" p={p:>3d}: F={F:.4f}, p={pval:.6f} {sig}")
# ─────────────────────────────────────────────────────────
# 7. PARTIAL CORRELATION
# ─────────────────────────────────────────────────────────
print(f"\n--- Partial Correlation (AR residuals) ---")
for ar_p in [12, 24, 60, 132]:
n = len(temp_series)
if ar_p >= n // 2: continue
try:
y_t = []; X_t = []
for t in range(ar_p, n):
y_t.append(temp_series[t])
X_t.append([1.0] + [temp_series[t-i] for i in range(1, ar_p+1)])
beta_t, _ = ols_fit(X_t, y_t)
resid_t = [y_t[i] - sum(X_t[i][j]*beta_t[j] for j in range(len(beta_t))) for i in range(len(y_t))]
y_s = []; X_s = []
for t in range(ar_p, n):
y_s.append(ssn_series[t])
X_s.append([1.0] + [ssn_series[t-i] for i in range(1, ar_p+1)])
beta_s, _ = ols_fit(X_s, y_s)
resid_s = [y_s[i] - sum(X_s[i][j]*beta_s[j] for j in range(len(beta_s))) for i in range(len(y_s))]
r = pearson_r(resid_t, resid_s)
n_eff = len(resid_t)
t_stat = r * math.sqrt((n_eff-2)/(1-r**2)) if abs(r) < 1 else float('inf')
print(f" AR({ar_p:>3d}): r={r:.4f}, t={t_stat:.3f}, n={n_eff}")
except: continue
# ─────────────────────────────────────────────────────────
# 8. PERMUTATION TEST
# ─────────────────────────────────────────────────────────
print(f"\n--- Permutation Test (p=12, 500 shuffles) ---")
def granger_f_only(y_series, x_series, p):
n = len(y_series)
if p >= n // 2: return None
y_dep = []; X_r = []; X_u = []
for t in range(p, n):
y_dep.append(y_series[t])
row_r = [1.0] + [y_series[t-i] for i in range(1, p+1)]
X_r.append(row_r)
row_u = row_r[:] + [x_series[t-j] for j in range(1, p+1)]
X_u.append(row_u)
try:
_, rss_r = ols_fit(X_r, y_dep)
_, rss_u = ols_fit(X_u, y_dep)
except: return None
dof_den = len(y_dep) - 2*p - 1
if dof_den <= 0 or rss_u <= 0: return None
return ((rss_r - rss_u) / p) / (rss_u / dof_den)
F_obs = granger_f_only(temp_series, ssn_series, 12)
count = 0; ssn_sh = ssn_series[:]
for _ in range(500):
random.shuffle(ssn_sh)
F_p = granger_f_only(temp_series, ssn_sh, 12)
if F_p is not None and F_p >= F_obs: count += 1
print(f" F_obs={F_obs:.4f}, perm_p={(count+1)/501:.4f}")
# ─────────────────────────────────────────────────────────
# 9. SLIDING WINDOW
# ─────────────────────────────────────────────────────────
print(f"\n--- Sliding Window (50yr, step=5yr, p=12) ---")
start = 0
while start + 600 <= N:
end = start + 600
res = granger_test(temp_series[start:end], ssn_series[start:end], 12)
if res:
F, pval, _, _, _ = res
yr_s = common_keys[start][0]; yr_e = common_keys[end-1][0]
sig = " *" if pval < 0.05 else ""
print(f" {yr_s}-{yr_e}: F={F:.4f}, p={pval:.6f}{sig}")
start += 60
# ─────────────────────────────────────────────────────────
# 10. DETRENDED ANALYSIS
# ─────────────────────────────────────────────────────────
print(f"\n--- Detrended Granger ---")
X_tr = [[1.0, i] for i in range(N)]
beta_tr, _ = ols_fit(X_tr, temp_series)
temp_dt = [temp_series[i] - (beta_tr[0] + beta_tr[1]*i) for i in range(N)]
print(f" Trend: {beta_tr[1]*120:.4f} °C/decade")
for p in all_lags:
res = granger_test(temp_dt, ssn_series, p)
if res:
F, pval, _, _, _ = res
sig = " *" if pval < 0.05 else ""
print(f" p={p:>3d}: F={F:.4f}, p={pval:.6f}{sig}")
# ─────────────────────────────────────────────────────────
# 11. NONLINEAR EXTENSIONS
# ─────────────────────────────────────────────────────────
print(f"\n--- Nonlinear: SSN² → Temp ---")
ssn_sq = [x**2 for x in ssn_series]
for p in [1, 6, 12, 24]:
res = granger_test(temp_series, ssn_sq, p)
if res:
F, pval, _, _, _ = res
sig = " *" if pval < 0.05 else ""
print(f" p={p:>3d}: F={F:.4f}, p={pval:.6f}{sig}")
def detrend(s):
n = len(s)
X = [[1.0, i] for i in range(n)]
beta, _ = ols_fit(X, s)
return [s[i] - (beta[0] + beta[1]*i) for i in range(n)]
mid = N // 2
print(f"\n--- Epoch-split (at {common_keys[mid][0]}) ---")
for label, t_sub, s_sub in [("First", temp_series[:mid], ssn_series[:mid]),
("Second", temp_series[mid:], ssn_series[mid:])]:
t_dt = detrend(t_sub)
for p in [6, 12, 24]:
res = granger_test(t_dt, s_sub, p)
if res:
F, pval, _, _, _ = res
sig = " *" if pval < 0.05 else ""
print(f" {label} half p={p:>3d}: F={F:.4f}, p={pval:.6f}{sig}")
# Regime-separated test
def rolling_mean(s, w=132):
return [sum(s[max(0,i-w+1):i+1]) / (i - max(0,i-w+1) + 1) for i in range(len(s))]
ssn_rm = rolling_mean(ssn_series, 132)
med_rm = sorted(ssn_rm)[len(ssn_rm)//2]
def contiguous_runs(indices, min_len=200):
runs = []; start = end = indices[0]
for i in range(1, len(indices)):
if indices[i] == end + 1: end = indices[i]
else:
if end - start + 1 >= min_len: runs.append((start, end))
start = end = indices[i]
if end - start + 1 >= min_len: runs.append((start, end))
return runs
high_idx = [i for i in range(N) if ssn_rm[i] >= med_rm]
low_idx = [i for i in range(N) if ssn_rm[i] < med_rm]
high_runs = contiguous_runs(high_idx) if high_idx else []
low_runs = contiguous_runs(low_idx) if low_idx else []
print(f"\n--- Regime-separated Granger ---")
for label, runs in [("HIGH", high_runs), ("LOW", low_runs)]:
for rs, re in runs:
t_run = detrend(temp_series[rs:re+1])
s_run = ssn_series[rs:re+1]
yr_s = common_keys[rs][0]; yr_e = common_keys[re][0]
n_run = re - rs + 1
for p in [6, 12, 24]:
if p >= n_run // 3: continue
res = granger_test(t_run, s_run, p)
if res:
F, pval, _, _, _ = res
sig = " *" if pval < 0.05 else ""
print(f" {label} {yr_s}-{yr_e} p={p:>3d}: F={F:.4f}, p={pval:.6f}{sig}")
# Absolute-residual correlation
print(f"\n--- Absolute-residual correlation ---")
for ar_p in [12, 24]:
n = len(temp_series)
try:
y_t = []; X_t = []
for t in range(ar_p, n):
y_t.append(temp_series[t])
X_t.append([1.0] + [temp_series[t-i] for i in range(1, ar_p+1)])
beta_t, _ = ols_fit(X_t, y_t)
resid_t = [abs(y_t[i] - sum(X_t[i][j]*beta_t[j] for j in range(len(beta_t)))) for i in range(len(y_t))]
y_s = []; X_s = []
for t in range(ar_p, n):
y_s.append(ssn_series[t])
X_s.append([1.0] + [ssn_series[t-i] for i in range(1, ar_p+1)])
beta_s, _ = ols_fit(X_s, y_s)
resid_s = [abs(y_s[i] - sum(X_s[i][j]*beta_s[j] for j in range(len(beta_s)))) for i in range(len(y_s))]
r = pearson_r(resid_t, resid_s)
print(f" AR({ar_p}): r(|resid_t|, |resid_s|) = {r:.4f}")
except: continue
print("\nAnalysis complete.")
```
## Expected Output Summary
- **N = 2,100** monthly observations (overlapping period of SILSO and NOAA NCEI records)
- **ADF tests**: SSN stationary (ADF = −8.20, p < 0.01); Temp non-stationary (ADF = −1.69); ΔTemp stationary (ADF = −19.86, p < 0.01)
- **Lag selection**: AIC = 8, BIC = 6 for temperature AR(p)
- **Granger (levels)**: No significance at any lag (smallest p = 0.169 at p=24)
- **Granger (differenced)**: No significance at any lag (smallest p = 0.127 at p=5)
- **Bidirectional**: No significance except marginal at p=6 (p=0.048, does not survive Bonferroni)
- **Partial correlations**: All |r| < 0.013
- **Permutation test** (p=12): p-value = 0.513
- **Sliding windows**: Only 1 of 26 significant (consistent with chance)
- **Detrended**: All null
- **Nonlinear (SSN²)**: One isolated nominal significance at p=6 (p=0.029), does not survive correction
- **Epoch-split**: All null in both halves
- **Regime-separated**: All null in all regimes
## Key Methods Implemented
1. **Augmented Dickey-Fuller test** with AIC-selected augmentation lags
2. **AIC/BIC lag selection** for AR model order
3. **Granger causality** on level-form and first-differenced series
4. **OLS via normal equations** with Gauss-Jordan matrix inversion (partial pivoting)
5. **F-distribution p-values** via regularized incomplete beta function
6. **Nonlinear extensions**: squared-SSN, epoch-split, regime-separated Granger tests
7. **Absolute-residual correlation** for variance-level coupling
8. **Permutation inference** as non-parametric validation
9. **Sliding-window temporal stability** analysis
10. **Bidirectional falsification** test
Discussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.