BEWARE FitAR requires R < 4.2 !!

We load libraries and the SOI datasets which is the best and longest intrumental proxy of ENSO.

suppressMessages(library(tseries))
suppressMessages(library(forecast))
suppressMessages(library(FitAR))
suppressMessages(library(strucchange)) 
vcov.ri <- function(x, ...) kernHAC(x, 
kernel = "Quadratic Spectral", prewhite = 1, approx = "AR(1)", ...)
#setwd("/Users/Shared/CCleaner")
raw <- read.csv("soi.csv")
summary(raw)
##       Date           SOI       
##  Min.   :1866   Min.   :-3850  
##  1st Qu.:1905   1st Qu.: -789  
##  Median :1945   Median : -126  
##  Mean   :1945   Mean   : -116  
##  3rd Qu.:1984   3rd Qu.:  551  
##  Max.   :2024   Max.   : 3314
DAT <- ts(raw, start=c(1866), end=c(2023),frequency=12)

The Chow test result clearly goes agaisnt the presence of a break (p-value 0.73).

sctest(raw$SOI ~ raw$Date, type = "Chow", point = 10)
## 
##  Chow test
## 
## data:  raw$SOI ~ raw$Date
## F = 0.31037, p-value = 0.7332

The ARMA model selection based on the AIC criterion suggest 14 lags. We test for a unit root with Phillips-Perron. As we obtain a very low p-value, we may thus reject the null hypothesis that there is one unit root, meaning that SOI is trend stationary (TS) over the entire period.

X = DAT[,"SOI"]
plot(X)

Acf(X)

SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)
##    p AIC-Exact AIC-Approx
## 1 14  25011.02  -1051.458
## 2 13  25018.67  -1050.759
## 3  4  25023.64  -1036.469
## 4  5  25025.10  -1034.469
## 5 11  25025.13  -1034.567
pp.test(X, alternative = c("stationary"))
## Warning in pp.test(X, alternative = c("stationary")): p-value smaller than
## printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  X
## Dickey-Fuller Z(alpha) = -856.07, Truncation lag parameter = 8, p-value
## = 0.01
## alternative hypothesis: stationary

Per the rpevious finding, we form the 14-months moving average.

Y1 = rollmean(DAT,14,align="center")
Y = Y1[,"SOI"]
plot(Y)

We test for a break in the Y trend using the first difference and applying a test for structural change. The Bayesian Information Criteria recommended by Bai & Perron suggests 2 breaks in September 1997 and December 1973 but with confidence intervals so large as to signify that no meaningful break exists.

U = diff(Y)
Z <- breakpoints(U ~ 1)
summary(Z) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = U ~ 1)
## 
## Breakpoints at observation number:
##                              
## m = 1                    1574
## m = 2               1289 1574
## m = 3           894 1289 1574
## m = 4       614 894 1289 1574
## m = 5   301 614 894 1289 1574
## 
## Corresponding to breakdates:
##                                                 
## m = 1                                    1997(9)
## m = 2                           1973(12) 1997(9)
## m = 3                   1941(1) 1973(12) 1997(9)
## m = 4           1917(9) 1941(1) 1973(12) 1997(9)
## m = 5   1891(8) 1917(9) 1941(1) 1973(12) 1997(9)
## 
## Fit:
##                                                          
## m   0        1        2        3        4        5       
## RSS 20838977 20810066 20779395 20758283 20713488 20704087
## BIC    22759    22771    22784    22797    22808    22822
plot(Z)

confid <- confint(Z,breaks = 2)
plot(U)
lines(Z)
lines(confid)
## Warning: Confidence intervals outside data time interval
##   from 1866(8) to 2022(6) (1871 observations)
## Warning: Overlapping confidence intervals

coef(Z, breaks = 2)
##                    (Intercept)
## 1866(8) - 1973(12)   -1.046769
## 1974(1) - 1997(9)    10.416792
## 1997(10) - 2022(6)   -9.727994
sapply(vcov(Z, breaks = 2, vcov = vcov.ri), sqrt)
## 1866(8) - 1973(12)  1974(1) - 1997(9) 1997(10) - 2022(6) 
##           5.878966          13.959590           9.901843
confint(Z, breaks = 2, vcov = vcov.ri)
## 
##   Confidence intervals for breakpoints
##   of optimal 3-segment partition: 
## 
## Call:
## confint.breakpointsfull(object = Z, breaks = 2, vcov. = vcov.ri)
## 
## Breakpoints at observation number:
##   2.5 % breakpoints 97.5 %
## 1 -3476        1289   4926
## 2   861        1574   3161
## 
## Corresponding to breakdates:
## Warning: Confidence intervals outside data time interval
##   from 1866(8) to 2022(6) (1871 observations)

## Warning: Overlapping confidence intervals
##   2.5 %    breakpoints 97.5 %  
## 1 1576(11) 1973(12)    2277(1) 
## 2 1938(4)  1997(9)     2129(12)
fac.ri <- breakfactor(Z, breaks = 1)
fm.ri <- lm(U ~ 0 + fac.ri)
plot(U)
lines(fitted(Z, breaks = 2), col = 2)

summary(fm.ri)
## 
## Call:
## lm(formula = U ~ 0 + fac.ri)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -339.10  -69.79   -0.46   69.49  328.54 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## fac.risegment1    1.029      2.660   0.387    0.699
## fac.risegment2   -9.728      6.123  -1.589    0.112
## 
## Residual standard error: 105.5 on 1869 degrees of freedom
## Multiple R-squared:  0.001429,   Adjusted R-squared:  0.0003601 
## F-statistic: 1.337 on 2 and 1869 DF,  p-value: 0.2629