Southern Oscillation Index Data

Description

The Southern Oscillation Index (SOI) is the difference in barometric pressure at sea level between Tahiti and Darwin. Annual SOI and Australian rainfall data, for the years 1900-2001, are given. Australia's annual mean rainfall is an area-weighted average of the total annual precipitation at approximately 370 rainfall stations around the country.

Usage

bomsoi

Format

This data frame contains the following columns:

Year

a numeric vector

Jan

average January SOI values for each year

Feb

average February SOI values for each year

Mar

average March SOI values for each year

Apr

average April SOI values for each year

May

average May SOI values for each year

Jun

average June SOI values for each year

Jul

average July SOI values for each year

Aug

average August SOI values for each year

Sep

average September SOI values for each year

Oct

average October SOI values for each year

Nov

average November SOI values for each year

Dec

average December SOI values for each year

SOI

a numeric vector consisting of average annual SOI values

avrain

a numeric vector consisting of a weighted average annual rainfall at a large number of Australian sites

NTrain

Northern Territory rain

northRain

north rain

seRain

southeast rain

eastRain

east rain

southRain

south rain

swRain

southwest rain

Source

Australian Bureau of Meteorology web pages:

http://www.bom.gov.au/climate/change/rain02.txt and http://www.bom.gov.au/climate/current/soihtm1.shtml

References

Nicholls, N., Lavery, B., Frederiksen, C. and Drosdowsky, W. 1996. Recent apparent changes in relationships between the El Nino – southern oscillation and Australian rainfall and temperature. Geophysical Research Letters 23: 3357-3360.

Examples

 
plot(ts(bomsoi[, 15:14], start=1900),
     panel=function(y,...)panel.smooth(1900:2005, y,...))
pause()

# Check for skewness by comparing the normal probability plots for 
# different a, e.g.
par(mfrow = c(2,3))
for (a in c(50, 100, 150, 200, 250, 300))
qqnorm(log(bomsoi[, "avrain"] - a))
  # a = 250 leads to a nearly linear plot

pause()

par(mfrow = c(1,1))
plot(bomsoi$SOI, log(bomsoi$avrain - 250), xlab = "SOI",
     ylab = "log(avrain = 250)")
lines(lowess(bomsoi$SOI)$y, lowess(log(bomsoi$avrain - 250))$y, lwd=2)
  # NB: separate lowess fits against time
lines(lowess(bomsoi$SOI, log(bomsoi$avrain - 250)))
pause()

xbomsoi <-
  with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
rainpos <- pretty(bomsoi$avrain, 5)
with(xbomsoi,
     {plot(cuberootRain ~ SOI, xlab = "SOI",
           ylab = "Rainfall (cube root scale)", yaxt="n")
     axis(2, at = rainpos^0.33, labels=paste(rainpos))
## Relative changes in the two trend curves
     lines(lowess(cuberootRain ~ SOI))
     lines(lowess(trendRain ~ trendSOI), lwd=2)
  })
pause()

xbomsoi$detrendRain <-
  with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
xbomsoi$detrendSOI <-
  with(xbomsoi, SOI - trendSOI + mean(trendSOI))
oldpar <- par(mfrow=c(1,2), pty="s")
plot(cuberootRain ~ SOI, data = xbomsoi,
     ylab = "Rainfall (cube root scale)", yaxt="n")
axis(2, at = rainpos^0.33, labels=paste(rainpos))
with(xbomsoi, lines(lowess(cuberootRain ~ SOI)))
plot(detrendRain ~ detrendSOI, data = xbomsoi,
  xlab="Detrended SOI", ylab = "Detrended rainfall", yaxt="n")
axis(2, at = rainpos^0.33, labels=paste(rainpos))
with(xbomsoi, lines(lowess(detrendRain ~ detrendSOI)))
pause()

par(oldpar)
attach(xbomsoi)
xbomsoi.ma0 <- arima(detrendRain, xreg=detrendSOI, order=c(0,0,0))
# ordinary regression model

xbomsoi.ma12 <- arima(detrendRain, xreg=detrendSOI,
                      order=c(0,0,12))
# regression with MA(12) errors -- all 12 MA parameters are estimated
xbomsoi.ma12
pause()

xbomsoi.ma12s <- arima(detrendRain, xreg=detrendSOI,
                      seasonal=list(order=c(0,0,1), period=12))
# regression with seasonal MA(1) (lag 12) errors -- only 1 MA parameter
# is estimated
xbomsoi.ma12s
pause()

xbomsoi.maSel <- arima(x = detrendRain, order = c(0, 0, 12),
                        xreg = detrendSOI, fixed = c(0, 0, 0,
                        NA, rep(0, 4), NA, 0, NA, NA, NA, NA),
                        transform.pars=FALSE)
# error term is MA(12) with fixed 0's at lags 1, 2, 3, 5, 6, 7, 8, 10
# NA's are used to designate coefficients that still need to be estimated
# transform.pars is set to FALSE, so that MA coefficients are not
# transformed (see help(arima))

detach(xbomsoi)
pause()

Box.test(resid(lm(detrendRain ~ detrendSOI, data = xbomsoi)),
          type="Ljung-Box", lag=20)

pause()

attach(xbomsoi)
 xbomsoi2.maSel <- arima(x = detrendRain, order = c(0, 0, 12),
                         xreg = poly(detrendSOI,2), fixed = c(0,
                         0, 0, NA, rep(0, 4), NA, 0, rep(NA,5)),
                         transform.pars=FALSE)
 xbomsoi2.maSel
qqnorm(resid(xbomsoi.maSel, type="normalized"))
detach(xbomsoi)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(DAAG)
Loading required package: lattice
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/DAAG/bomsoi.Rd_%03d_medium.png", width=480, height=480)
> ### Name: bomsoi
> ### Title: Southern Oscillation Index Data
> ### Aliases: bomsoi
> ### Keywords: datasets
> 
> ### ** Examples
>  
> plot(ts(bomsoi[, 15:14], start=1900),
+      panel=function(y,...)panel.smooth(1900:2005, y,...))
> pause()
> 
> # Check for skewness by comparing the normal probability plots for 
> # different a, e.g.
> par(mfrow = c(2,3))
> for (a in c(50, 100, 150, 200, 250, 300))
+ qqnorm(log(bomsoi[, "avrain"] - a))
>   # a = 250 leads to a nearly linear plot
> 
> pause()
> 
> par(mfrow = c(1,1))
> plot(bomsoi$SOI, log(bomsoi$avrain - 250), xlab = "SOI",
+      ylab = "log(avrain = 250)")
> lines(lowess(bomsoi$SOI)$y, lowess(log(bomsoi$avrain - 250))$y, lwd=2)
>   # NB: separate lowess fits against time
> lines(lowess(bomsoi$SOI, log(bomsoi$avrain - 250)))
> pause()
> 
> xbomsoi <-
+   with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
> xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
> xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
> rainpos <- pretty(bomsoi$avrain, 5)
> with(xbomsoi,
+      {plot(cuberootRain ~ SOI, xlab = "SOI",
+            ylab = "Rainfall (cube root scale)", yaxt="n")
+      axis(2, at = rainpos^0.33, labels=paste(rainpos))
+ ## Relative changes in the two trend curves
+      lines(lowess(cuberootRain ~ SOI))
+      lines(lowess(trendRain ~ trendSOI), lwd=2)
+   })
> pause()
> 
> xbomsoi$detrendRain <-
+   with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
> xbomsoi$detrendSOI <-
+   with(xbomsoi, SOI - trendSOI + mean(trendSOI))
> oldpar <- par(mfrow=c(1,2), pty="s")
> plot(cuberootRain ~ SOI, data = xbomsoi,
+      ylab = "Rainfall (cube root scale)", yaxt="n")
> axis(2, at = rainpos^0.33, labels=paste(rainpos))
> with(xbomsoi, lines(lowess(cuberootRain ~ SOI)))
> plot(detrendRain ~ detrendSOI, data = xbomsoi,
+   xlab="Detrended SOI", ylab = "Detrended rainfall", yaxt="n")
> axis(2, at = rainpos^0.33, labels=paste(rainpos))
> with(xbomsoi, lines(lowess(detrendRain ~ detrendSOI)))
> pause()
> 
> par(oldpar)
> attach(xbomsoi)
> xbomsoi.ma0 <- arima(detrendRain, xreg=detrendSOI, order=c(0,0,0))
> # ordinary regression model
> 
> xbomsoi.ma12 <- arima(detrendRain, xreg=detrendSOI,
+                       order=c(0,0,12))
> # regression with MA(12) errors -- all 12 MA parameters are estimated
> xbomsoi.ma12

Call:
arima(x = detrendRain, order = c(0, 0, 12), xreg = detrendSOI)

Coefficients:
         ma1      ma2     ma3      ma4      ma5     ma6      ma7     ma8
      0.0227  -0.1371  -0.031  -0.2511  -0.0710  0.0361  -0.0798  0.0559
s.e.  0.1027   0.1062   0.110   0.1193   0.1127  0.1413   0.1208  0.1226
         ma9     ma10     ma11     ma12  intercept  detrendSOI
      0.2163  -0.0075  -0.2933  -0.4603     7.5242      0.0353
s.e.  0.1666   0.1149   0.1563   0.1262     0.0068      0.0044

sigma^2 estimated as 0.08423:  log likelihood = -22.66,  aic = 75.32
> pause()
> 
> xbomsoi.ma12s <- arima(detrendRain, xreg=detrendSOI,
+                       seasonal=list(order=c(0,0,1), period=12))
> # regression with seasonal MA(1) (lag 12) errors -- only 1 MA parameter
> # is estimated
> xbomsoi.ma12s

Call:
arima(x = detrendRain, seasonal = list(order = c(0, 0, 1), period = 12), xreg = detrendSOI)

Coefficients:
         sma1  intercept  detrendSOI
      -0.4297     7.5223      0.0390
s.e.   0.0914     0.0191      0.0045

sigma^2 estimated as 0.09898:  log likelihood = -29.05,  aic = 66.11
> pause()
> 
> xbomsoi.maSel <- arima(x = detrendRain, order = c(0, 0, 12),
+                         xreg = detrendSOI, fixed = c(0, 0, 0,
+                         NA, rep(0, 4), NA, 0, NA, NA, NA, NA),
+                         transform.pars=FALSE)
> # error term is MA(12) with fixed 0's at lags 1, 2, 3, 5, 6, 7, 8, 10
> # NA's are used to designate coefficients that still need to be estimated
> # transform.pars is set to FALSE, so that MA coefficients are not
> # transformed (see help(arima))
> 
> detach(xbomsoi)
> pause()
> 
> Box.test(resid(lm(detrendRain ~ detrendSOI, data = xbomsoi)),
+           type="Ljung-Box", lag=20)

	Box-Ljung test

data:  resid(lm(detrendRain ~ detrendSOI, data = xbomsoi))
X-squared = 34.455, df = 20, p-value = 0.02321

> 
> pause()
> 
> attach(xbomsoi)
>  xbomsoi2.maSel <- arima(x = detrendRain, order = c(0, 0, 12),
+                          xreg = poly(detrendSOI,2), fixed = c(0,
+                          0, 0, NA, rep(0, 4), NA, 0, rep(NA,5)),
+                          transform.pars=FALSE)
>  xbomsoi2.maSel

Call:
arima(x = detrendRain, order = c(0, 0, 12), xreg = poly(detrendSOI, 2), transform.pars = FALSE, 
    fixed = c(0, 0, 0, NA, rep(0, 4), NA, 0, rep(NA, 5)))

Coefficients:
      ma1  ma2  ma3      ma4  ma5  ma6  ma7  ma8     ma9  ma10     ma11
        0    0    0  -0.4350    0    0    0    0  0.3888     0  -0.2950
s.e.    0    0    0   0.1089    0    0    0    0  0.1317     0   0.1174
         ma12  intercept       1       2
      -0.6759     7.5216  2.3964  0.2173
s.e.   0.1151     0.0068  0.2310  0.2480

sigma^2 estimated as 0.06832:  log likelihood = -23.18,  aic = 62.36
> qqnorm(resid(xbomsoi.maSel, type="normalized"))
> detach(xbomsoi)
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>