Time Series Analysis - Assignment 2

Time Series Analysis of Antarctica land ice mass series

Setup

# Clear the environment and the console
rm(list = ls()); cat("\f")



#----------------------------#
# Install required libraries #
#----------------------------#
packages <- c("ggplot2", "ggpubr", "readxl", "readr", "dplyr", "tidyr", "psych", "stringr",
              "lubridate", "knitr", "outliers", "MVN", "TSA", "tseries", "lmtest", "FSAdata",
              "forecast", "matrixcalc", "car", "corpcor", "scales", "QuantPsyc",
              "urca", "rugarch", "fGarch", "tswge",
              "imputeTS", "shiny", "coda", "rjags", "runjags", "ks", "epiDisplay", "fastDummies")
# Install any packages not already installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) { install.packages(packages[!installed_packages]) }
#-------------------------#
# Load libraries          #
#-------------------------#
invisible(lapply(packages, library, character.only = TRUE))
#-------------------------#
# Clear the environment and the console
rm(list = ls()); cat("\f")

Required packages

library(TSA)
library(readr)
library(stringr)
library(dplyr)
library(lubridate)

Introduction

Aims/Objectives

This report summarises the analysis of the Antarctica land ice mass series, and the process followed to identify suitable ARIMA(p,d,q) models as part of Assignment 2 for Time Series Analysis. The report covers analysis of the dataset, addressing non-stationarity, and using model specification tools to suggest possible ARIMA models.

Methodology

Creating a function for producing ACF and PACF plots

To reduce repetition of R codes, a function was created to produce ACF and PACF plots.

acf_pacf <- function(ts_object, plot_note){
 # Produce ACF and PACF plots for a time series object.
 # Optional: specify a plot note (as a string) to improve plot titles.
 # Example: acf_pacf(price_ts, "Box-Cox transformed, first differenced price series.")
 par(mfrow=c(1,2))
 if (missing(plot_note)) {
 acf(ts_object, main ="ACF plot", lag.max = min(c(length(ts_object), 50)))
 pacf(ts_object, main ="PACF plot", lag.max = min(c(length(ts_object), 50)))
 } else {
 acf(ts_object, main=paste("ACF plot of", plot_note), lag.max = min(c(length(ts_object),50)))
 pacf(ts_object, main=paste("PACF plot of", plot_note), lag.max = min(c(length(ts_object),50)))
 }
 par(mfrow=c(1,1))
}

Retrieving the data

A summary of the Antarctica land ice mass data imported into R showed the imported dataset had 19 observations, including a minimum date of 2002 and a maximum date of 2020, the same as the original data file.

ice <- read_csv("assignment2Data2022.csv", col_names = TRUE)
colnames(ice) <- c("Year", "Annual_ice_mass")
describe(ice, fast=TRUE, quant=c(.25, .5, .75))

The data summary also showed that between 2002 and 2020, the mean annual change in land ice mass relative to 2001 was -1,063.77 billion metric tons. Considering the timeframe for the series was 19 years, it was shocking to note that the ‘Annual_ice_mass’ had a range of 2,531.51 billion metric tons, more than double the mean annual change and indicative of a very large change in land ice mass over the course of less than 20 years.

As the Antarctic land Ice dataset had been stored as a dataframe when it was first imported into R, it was converted to a time series object with a frequency of 1 as each observation was an annual change.

ice_ts <- ts(ice$Annual_ice_mass, start = 2002, frequency = 1)
ice_ts

Time Series:
Start = 2002 
End = 2020 
Frequency = 1 
 [1]    -7.546667  -115.186364  -262.302500  -228.376667  -129.987500  -316.476667
 [7]  -587.155833  -562.316667  -840.778333  -938.991111 -1073.690909 -1283.981111
[13] -1494.724444 -1846.245556 -1768.321111 -1801.952000 -2146.672000 -2267.945000
[19] -2539.055000

Exploring the data

The Antarctic land ice series was plotted to visually inspect the series for the 5 key time series traits; trend, seasonality, changing variance, behaviour, and change point/intervention.

plot(ice_ts, type="o", xlab="Year",
 ylab="Change in Antarctica land ice mass, relative to 2001 (billion metric tons)",
 main="Time series plot of the NASA Annual Antarctica land ice mass series.")

Figure 1: Time series plot of the Annual Antarctica land ice mass series.

The time series plot of the Antarctica land ice mass series showed the following time series traits;

1. Trend

From the time series plot in Figure 1, it was clear the series had a downward, linear trend, indicating that the Antarctic land ice mass was getting smaller relative to 2001 as time progressed.

2. Seasonality

There were no clear signs of seasonality in Figure 1.

3. Changing variance

There were no clear signs of changing variance in Figure 1.

4. Behaviour

The time series plot in Figure 1 exhibited autoregressive (AR) behaviour, with many successive points. There were no clear fluctuations in the series, so there was no evidence of moving average (MA) behaviour.

5. Change point

There was no clear evidence of a change point (intervention) in the Antarctica land ice mass series as the series appeared to exhibit the same trend and behaviour, with no sudden changes. This was not surprising considering the factors that would affect Antarctic land ice mass (e.g., climate change) do not change substantially from year to year. It would take a very large intervention (e.g., a meteor strike) to change the behaviour of the series.

The points in Figure 1 show evidence of succeeding measurements being related to one another, and Figure 2 highlights this relationship more clearly using a scatter plot of neighbouring pairs.

#lag the time series
par(mfrow=c(1,1))
plot(y=ice_ts, x=zlag(ice_ts), 
 ylab="Change in Antarctica land ice mass, relative to 2001 (billion metric tons)", 
 xlab="Change in Antarctica land ice mass in previous year",
 main = "Scatter plot of Change in Antarctica land ice mass in consecutive Years.")

Figure 2: Scatter plot of each Annual Antarctica land ice mass measurement and the measurement taken in the previous year.

# correlation of neighbouring pairs
y = ice_ts # Assign the ice data to y
x = zlag(ice_ts) # Generate first lag of the ice series
index = 2:length(x) # Create an index to get rid of the first NA value in x
cor(y[index],x[index]) # Calculate correlation between numerical values in x and y

[1] 0.985486

The scatter plot in Figure 2 shows a strong, upward trend, indicative of a strong correlation between neighbouring pairs. The plot indicates that small annual changes in Antarctic land ice mass (relative to 2001) tend to be followed by small changes, moderate changes tend to be followed by moderate changes, and large changes tend to be followed by large changes. These observations were supported by the correlation between neighbouring annual land ice mass changes which, at 0.985, indicated a very strong, positive correlation between neighbouring points.

Assessing stationarity

Figure 1 showed a strong, downward, linear trend in the Antarctica ice series, suggesting that the raw Antarctica ice series was non-stationary. As ARIMA models can only be produced using a stationary series, ACF and PACF plotswere produced to assess the series for stationarity.

acf_pacf(ice_ts, "the Antarctica Ice series.")

Figure 3: ACF and PACF plots of the Antarctica land ice mass series.

The ACF plot in Figure 3 showed 3 significant autocorrelation lags and a decaying pattern, indicative of a non-stationary, autoregressive process. There was no evidence of seasonality in the series as there was no clear wave pattern in the ACF lags.

The first partial autocorrelation lag in the PACF plot in Figure 3 was highly significant, while all other lags were not significant, indicative of a non-stationary, autoregressive series. Several unit root tests, including an Augmented Dickey-Fuller test, were also used to assess the series for stationarity, with results shown in Table 4.

adf.test(ice_ts, alternative = c("stationary")) # significant=stationary, doubt range:(0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  ice_ts
Dickey-Fuller = -2.7141, Lag order = 2, p-value = 0.3004
alternative hypothesis: stationary

pp.test(ice_ts) # significant = stationary

    Phillips-Perron Unit Root Test

data:  ice_ts
Dickey-Fuller Z(alpha) = -5.8754, Truncation lag parameter = 2, p-value = 0.7516
alternative hypothesis: stationary

kpss.test(ice_ts) # significant = non-stationary

    KPSS Test for Level Stationarity

data:  ice_ts
KPSS Level = 0.72818, Truncation lag parameter = 2, p-value = 0.01098

Table 4: Statistical tests for stationarity and normality in the Antarctica land ice mass series.

The null hypothesis under the Augmented Dickey-Fuller test and the Phillips-Perron (PP) test is that the series Is non-stationary. With p-values of 0.30 for the Augmented Dickey-Fuller test, and 0.75 for the PP test, there was insufficient evidence to reject the null hypothesis of both tests and therefore we could conclude that the series was not stationary. This conclusion was supported by results of a KPSS test, which produced a p-value of 0.01, indicating there was sufficient evidence to reject the null hypothesis under the KPSS test that the series was stationary. Finally, a Q-Q plot was produced to assess the raw series for normality, alongside a Shapiro-Wilk test. Although the Q-Q plot in Figure 4 showed deviations from the diagonal in each tail of the series, the Shapiro-Wilk test produced a p-value of 0.147, indicating that the measurements in the Antarctica land ice mass series were approximately normally distributed.

## Normality checks
qqnorm(ice_ts)
qqline(ice_ts, col = 2)

Figure 4: Normal Q-Q plot of the Antarctica land ice mass series.

shapiro.test(ice_ts) # not significant - approximately normal

    Shapiro-Wilk normality test

data:  ice_ts
W = 0.92616, p-value = 0.1471

Table 5: Result of the Shapiro-Wilk test on the raw Antarctica Ice series.

Addressing non-stationarity

Transformation

Although there was no clear evidence of changing variance in Figure 1, a Box-Cox transformation was applied to the series to observe the impact of transformation on the series.

As the series contained negative values (declines in the land ice mass relative to 2001), a constant had to be added to the series to make all values greater than zero before a Box-Cox transformation could be applied. A vector of candidate power transformation values was also required to obtain candidate lambda values for the Box-Cox transformation due to computational limits. Through trial and error, limits of zero and 2 were identified as appropriate power transformation limits for the Box-Cox transformation.

The plot of candidate lambda values for the Box-Cox transformation of the Antarctica land ice mass series can be seen in Figure 5. Table 6 shows the values of the first and third vertical lines from Figure 5 were 0.65 and 1.22, respectively, and the chosen lambda value for the Box-Cox transformation (the middle vertical line in Figure 5) was 0.89.

#BC = BoxCox.ar(ice_ts) # ERROR: data values must be positive.
min(ice_ts) # minimum is -2539.055

[1] -2539.055

# add the minimum to each value in the series, to make all values greater than zero.
ice_ts2 <- ice_ts + abs(min(ice_ts)) + 0.01
BC <- BoxCox.ar(y=ice_ts2, lambda=seq(0, 2, 0.01))

Figure 5: Log likelihood plot for the lambda values in the Box-Cox transformation of the Antarctica land ice mass series.

BC$ci #Values of the first and third vertical lines

[1] 0.65 1.22

# To find the lambda value of the middle vertical line
lambda <- BC$lambda[which(max(BC$loglike) == BC$loglike)]
lambda

[1] 0.89

Table 6: Results for candidate lambda values from the Box-Cox power transformation function.

With an appropriate lambda value identified, the Box-Cox transformation was applied to the Antarctica Ice series and the resultant series was plotted for visual inspection in Figure 6.

# Apply Box-Cox transformation
iceBC <- ((ice_ts2^lambda) - 1) / lambda
plot(iceBC, type="o", xlab="Years", 
 ylab="Change in Antarctica land ice mass, relative to 2001 (billion metric tons)",
 main = "Time series plot of BoxCox transformed Antarctica Ice series.")

Figure 6: Time series plot of the Box-Cox transformed Annual Antarctica land ice mass series.

The time series plot in Figure 6 showed no clear change or improvement in the variation of the series, and the Q-Q plot in Figure 7 still showed deviations from the diagonal at both tails of the series. There was also little change in the result from the Shapiro-Wilk test of the Box-Cox transformed series, shown in Table 7, indicating that values in the series were still approximately normally distributed.

qqnorm(iceBC, main = "Normal Q-Q Plot for the Box-Cox transformed Antarctica Ice series.")
qqline(iceBC, col = 2)

Figure 7: Normal Q-Q Plot for the Box-Cox transformed Antarctica Ice series.

shapiro.test(iceBC)

    Shapiro-Wilk normality test

data:  iceBC
W = 0.92582, p-value = 0.145

Table 7: Results from the Shapiro-Wilk test of the Box-Cox transformed Antarctic Ice series.

Finally, the Augmented Dickey-Fuller test (Table 8) produced a p-value of 0.563, indicating that the Box-Cox transformed Antarctica Ice series was non-stationary.

adf.test(iceBC) # significant = stationary, doubt range: (0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  iceBC
Dickey-Fuller = -2.0237, Lag order = 2, p-value = 0.5634
alternative hypothesis: stationary

Table 8: Results from the Augmented Dickey-Fuller test of the Box-Cox transformed Antarctic Ice series.

As there was no evidence to suggest the raw Antarctica Ice series exhibited changing variation, and the Box-Cox transformation did not change or improve the normality of the series, the Box-Cox transformation was deemed unnecessary, and the raw series would be used when applying differencing to address non-stationarity.

Differencing

Differencing was applied to the raw Antarctica Ice series to address the trend in the raw series, with the resulting series plotted in Figure 8.

ice_diff <- diff(ice_ts, differences = 1)
par(mar=c(5,6,4,1)+.1) #set plot window margins
plot(ice_diff, type="o", xlab="Year",
 ylab="Change in Antarctica land ice mass, relative to 2001 \n(billion metric tons)",
 main="Time series plot of first differenced Antarctica Ice series.")

Figure 8: Time series plot of first differenced Antarctica Ice series

The time series plot in Figure 8 appeared to exhibit a slight downward trend, indicating that the first differenced series may still be non-stationary. However, the ACF and PACF plots of the first differenced Antarctica Ice series (Figure 9) showed no significant points and no decaying pattern, suggesting the first differenced series may have been stationary.

acf_pacf(ice_diff, "the first differenced \nAntarctica Ice series.")

Figure 9: ACF and PACF plots of the first differenced Antarctica land ice mass series.

To test for stationarity, an Augmented Dickey-Fuller test of the first differenced series was performed, with the test resulting in a p-value of 0.357 which, at the α=0.05 level, indicated that the first differenced Antarctica Ice series was non-stationary.

adf.test(ice_diff)# significant = stationary, doubt range: (0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  ice_diff
Dickey-Fuller = -2.5664, Lag order = 2, p-value = 0.3566
alternative hypothesis: stationary

Table 9: Result from the Augmented Dickey-Fuller test of the first differenced Antarctic Ice series.

As first differencing had not produced a stationary series, second differencing was applied to the raw Antarctica Ice series, with the resulting series plotted in Figure 10.

ice_diff2 <- diff(ice_ts, differences = 2)
plot(ice_diff2, type="o", xlab="Year",
 ylab="Change in Antarctica land ice mass, relative to 2001 \n(billion metric tons)",
 main="Time series plot of second differenced Antarctica Ice series.")

Figure 10: Time series plot of second differenced Antarctica Ice series.

The second differenced time series in Figure 10 showed no clear trend, suggesting that second differencing may have resulted in a stationary series.

acf_pacf(ice_diff2, "the \nsecond differenced Antarctica Ice series")

Figure 11: ACF and PACF plots of the second differenced Antarctica land ice mass series.

The ACF and PACF plots of the second differenced series (Figure 11) showed a significant PACF lag at the second lag, suggesting the second differenced series may be non-stationary. Supporting this was the Augmented Dickey-Fuller test of the second differenced series (Table 10), which produced a p-value of 0.073, indicating that, at the α=0.05 level, the second differenced series may still have been non-stationary.

However, as the result of the Augmented Dickey-Fuller test was within the “doubt range” of 0.03 to 0.1, additional unit root tests were performed to better evaluate whether second differencing had resulted in a stationary series.

adf.test(ice_diff2)# significant = stationary, doubt range: (0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  ice_diff2
Dickey-Fuller = -3.438, Lag order = 2, p-value = 0.0725
alternative hypothesis: stationary

Table 10: Results from the Augmented Dickey-Fuller test, Phillips-Perron test, and KPSS test of the second differenced Antarctic Ice series.

A Phillips-Perron Unit Root test and a KPSS test (Table 11) produced p-values of 0.026, and 0.1, respectively, indicating there was sufficient evidence to conclude the second differenced series was stationary.

# p-value is in the doubt range, so refer to additional tests.
pp.test(ice_diff2) # significant = stationary

    Phillips-Perron Unit Root Test

data:  ice_diff2
Dickey-Fuller Z(alpha) = -19.817, Truncation lag parameter = 2, p-value = 0.02604
alternative hypothesis: stationary

kpss.test(ice_diff2) # significant = non-stationary

    KPSS Test for Level Stationarity

data:  ice_diff2
KPSS Level = 0.11554, Truncation lag parameter = 2, p-value = 0.1

Table 11: Results from the Phillips-Perron test, and KPSS test of the second differenced Antarctic Ice series.

As the series appeared stationary in Figure 10, and the Phillips-Perron Unit Root test and KPSS test in Table 11 indicated the second differenced series was stationary, it was concluded that second differencing of the Antarctica Ice series had resulted in a stationary series.

As second differencing the Antarctica Ice series had resulted in a stationary series, it was deemed suitable for use in model specification. All proposed ARIMA(p,d,q) models would have a value of “d” equal to 2, as second differencing had been applied.

Model specification

ACF & PACF plots

The first step in identifying potential models for the Antarctica Ice series was to use ACF and PACF plots to propose values for the orders of autoregression (p) and moving average (q).

acf_pacf(ice_diff2, "the \nsecond differenced Antarctica Ice series")

Figure 12: ACF and PACF plots of the second differenced Antarctica Ice series.

From the PACF plot in Figure 12, the second partial autocorrelation lag was deemed significant, and the first lag was deemed borderline significant, indicating that the value of “p” could be 1 or 2. These values were understandable given autoregressive behaviour was identified in the raw Antarctica Ice series.

There were no clearly significant autocorrelation lags in the ACF plot of the second differenced series. However, the first lag of the ACF plot was deemed borderline significant, and thus the value of “q” was proposed as either 0 or 1.

The resulting proposed models from the ACF and PACF plots were:

• { ARIMA(1,2,0), ARIMA(2,2,0), ARIMA(1,2,1), ARIMA(2,2,1) }

EACF

The second step in identifying potential models for the Antarctica Ice series involved using the extended ACF function, or EACF. Running the EACF function on the second differenced series with default parameters resulted in an error, which was overcome by specifying the maximum values for AR and MA in the EACF function. The resulting output is shown in Table 12.

# eacf(ice_diff2) #ERROR: requires numeric/complex matrix/vector arguments
# eacf(ice_diff2, ar.max = 5, ma.max = 5) #ERROR: requires numeric/complex matrix/vector arguments
eacf(ice_diff2, ar.max = 4, ma.max = 3)

AR/MA
  0 1 2 3
0 o o o o
1 x x o o
2 o x o o
3 o o o o
4 o o o o

Table 12: EACF plot of the second differenced Antarctic Ice series.

The most top-left point in the EACF plot was selected as (0,0), and the models proposed using the EACF plot were:

• { ARIMA(0,2,0), ARIMA(0,2,1) }

It should be noted that the ARIMA(0,2,0) model will have no coefficients. Trend modelling would be used to see if this model captures the autocorrelation in the series.

BIC

The third step in identifying potential models for the Antarctica Ice series was to use the Bayesian Information Criterion, or “BIC”, to propose values for “p” and “q”.

The BIC plot was first produced using values of 10 for ‘nma’ and ‘nar’ however, this resulted in an error. As it was known from the ACF, PAC, and EACF plots that the orders of “p” and “q” would likely be small, the values of ‘nma’ and ‘nar’ were reduced to 5, and then to 3, with the resulting BIC plots shown in Figure 13 and Figure 14, respectively.

#res = armasubsets(y=ice_diff2, nar=10, nma=10, y.name='p', ar.method='ols') #ERRORs
#plot(res) #ERRORs
## evidence so far suggests p and q will be small, so try with a reduced value for 'nma' and 'nar'.
res2 = armasubsets(y=ice_diff2, nar=5, nma=5, y.name='p', ar.method='ols')
plot(res2) # algorithm hits a BIC value that was NaN

Figure 13: BIC plot of the second differenced Antarctic Ice series, using ‘nar’ and ‘nma’ values of 5.

# evidence so far suggests p and q will be small, so try with a further reduced value for 'nma' and 'nar'.
res3 = armasubsets(y=ice_diff2, nar=3, nma=3, y.name='p', ar.method='ols')
plot(res3) # algorithm hits a BIC value that was inf

Figure 14: BIC plot of the second differenced Antarctic Ice series, using ‘nar’ and ‘nma’ values of 3.

The BIC plot with values of 3 for ‘nar’ and ‘nma’ (Figure 14) proposed models that were in-line with the ACF, PACF, and EACF plots, so it was selected for use in proposing models for the second differenced series.

The model proposals from the BIC model specification method were:

• { ARIMA(1,2,1), ARIMA(2,2,1), ARIMA(3,2,1), ARIMA(1,2,0), ARIMA(2,2,0) }

Results

Final set of possible models

The final set of possible models for the Antarctica land ice mass series was:

• { ARIMA(1,2,0),
• ARIMA(2,2,0),
• ARIMA(1,2,1),
• ARIMA(2,2,1),
• ARIMA(0,2,0),
• ARIMA(0,2,1),
• ARIMA(3,2,1) }

Most of these models made sense in the context of the data exploration, which identified autoregressive behaviour in the raw Antarctica Ice series (i.e., “p” would likely be non-zero) and no clear evidence of moving average behaviour (i.e., “q” would likely be low or zero). It should be noted that the ARIMA(0,2,0) model will have no coefficients, and trend modelling would be required to see if this model captures the autocorrelation in the series.

The following models were proposed by more than one model specification method;

• { ARIMA(1,2,0),
• ARIMA(2,2,0),
• ARIMA(1,2,1),
• ARIMA(2,2,1) }

Further testing, including diagnostic checking, would be required to identify the optimal model from the set of proposed models.

Summary & Conclusion

The objective of the analysis and model specification conducted was to understand the Antarctica land ice mass dataset, and to propose a set of possible ARIMA(p,d,q) models for the series.

Exploration of the Antarctica land ice mass dataset revealed a large range and a mean change of -1,063.77 billion metric tons, relative to the ice mass in 2001.

Analysis of the Antarctica Ice series identified that the series exhibited a clear downward trend and autoregressive behaviour (successive points), with no clear signs of seasonality, changing variance, or a change point. A high correlation between neighbouring points was also identified, which aligned with the autocorrelation behaviour observed when the series was plotted.

ACF and PACF plots of the Antarctica Ice series identified significant autocorrelation lags, and a single, very significant partial autocorrelation lag, indicative of a non-stationary series. Further tests confirmed that the raw series was not stationary, but that the measurements in the series were approximately normally distributed.

Although the series did not exhibit changing variance, a Box-Cox transformation was applied to the series to observe its impact. Visual inspection, a Shapiro-Wilk test, and a unit root test identified that the Box-Cox transformation did not change the series much, so the transformation was deemed unnecessary.

Differencing was then applied to the raw series to address the trend in the series. Visual inspection and a unit root test confirmed that first differencing did not overcome the non-stationarity, and so second differencing was applied. Although the Augmented Dickey-Fuller test of the second differenced series suggested the series was still non-stationary, visual inspection, and further unit root tests (PP and KPSS) confirmed that second differencing had made the series stationary. As such, the second differenced series was chosen for use in model specification.

ACF, PACF, EACF and BIC model specification tools were used to proposed suitable ARIMA(p,d,q) models for the Antarctica Ice series, with the final set of proposed models containing 7 ARIMA(p,d,q) models;

{ ARIMA(1,2,0), ARIMA(2,2,0), ARIMA(1,2,1), ARIMA(2,2,1), ARIMA(0,2,0), ARIMA(0,2,1), ARIMA(3,2,1) }.