Introduction

The dataset use in this project describes the minimum daily temperature over 10 years (19981 and 1990) in Melbourne, Australia. The units are in degree celsius and there are 3650 observations. The source of the data was created by Australian Bureau of meteorology. We will analyze this dataset to capture the trend, applied transformation on the data where required, plot the ACF & PACF’s and coefficient tests on suitable models, and forecast for the next five years.

Require packages

# loading relevant packages into rstudio 
library(readr)
library(dplyr)
library(tidyr)
library(TSA)
library(tseries)
library(lmtest)
library(forecast)
library(fUnitRoots)
library(FitAR)
library(astsa)

Import Dataset Convert To Time Series

# reading the dataset  into R

Mel_Min_temp <- read.csv("Minimun Daily Temperature.csv", header = TRUE)


head(Mel_Min_temp,5)
##   ï..Year Month Day Minimum.Daily.Temperature.Degree..Degree.C..
## 1    1981     1   1                                         20.7
## 2    1981     1   2                                         17.9
## 3    1981     1   3                                         18.8
## 4    1981     1   4                                         14.6
## 5    1981     1   5                                         15.8
# converting the dataframe to time series object

min_TempTS <- ts(Mel_Min_temp$Minimum.Daily.Temperature.Degree..Degree.C..,
start=c(1981, 1), end=c(1990, 365), frequency = 365)

#checking the class of the dataset 
class(min_TempTS)
## [1] "ts"

1 Descriptive Analysis

#plotting time Series data 

plot(min_TempTS , col=c("blue"), ylab='Min temperature', xlab='Time', type='o',
main = "Time series plot of minimum temperature in Melbourne.")

Figure 1: shows the changes in daily minimum temperature in melbourne from 1981 to 1990. We can conclude the plot as well as some strong cyclic behaviour. there is no apparent trend in the data over this period.

#converting daily minimum temperature to a monthly mean teamperature to better see seasonality 

monthlyMinTemp <- aggregate(Mel_Min_temp$Minimum.Daily.Temperature.Degree..Degree.C..~Month+Mel_Min_temp$ï..Year,
Mel_Min_temp, mean )

head(monthlyMinTemp,5)
##   Month Mel_Min_temp$ï..Year
## 1     1                 1981
## 2     2                 1981
## 3     3                 1981
## 4     4                 1981
## 5     5                 1981
##   Mel_Min_temp$Minimum.Daily.Temperature.Degree..Degree.C..
## 1                                                 17.712903
## 2                                                 17.678571
## 3                                                 13.500000
## 4                                                 12.356667
## 5                                                  9.490323
# converting monthly temperature data  into times series  data

monthlyMinTempTS <- ts(monthlyMinTemp$`Mel_Min_temp$Minimum.Daily.Temperature.Degree..Degree.C..`,
start=c(1981, 1), end=c(1990, 12), frequency = 12)

monthlyMinTempTS
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1981 17.712903 17.678571 13.500000 12.356667  9.490323  7.306667  7.577419
## 1982 16.567742 15.921429 14.935484 11.470000  9.583871  5.606667  4.641935
## 1983 13.180645 16.807143 15.777419 10.596667 10.116129  6.600000  6.890323
## 1984 14.309677 14.944828 12.867742 10.750000  8.112903  7.730000  5.987097
## 1985 14.219355 14.032143 15.877419 12.976667  9.419355  7.073333  6.135484
## 1986 13.825806 14.196429 14.690323 11.653333 10.274194  7.526667  6.961290
## 1987 13.235484 13.889286 12.619355 12.250000  9.806452  8.273333  5.983871
## 1988 16.493548 14.524138 14.748387 12.833333 11.387097  8.386667  8.232258
## 1989 15.180645 16.371429 15.803226 12.563333 10.725806  6.560000  6.332258
## 1990 15.577419 15.417857 14.835484 13.433333  9.748387  7.720000  8.183871
##            Aug       Sep       Oct       Nov       Dec
## 1981  7.238710 10.143333 10.087097 11.890000 13.680645
## 1982  7.903226  7.280000  9.545161 12.486667 13.754839
## 1983  8.706452  9.210000 10.312903 11.993333 14.396774
## 1984  8.696774  8.046667 10.632258 12.623333 12.643333
## 1985  7.635484  8.803333 10.490323 13.073333 14.109677
## 1986  7.387097  8.933333  9.683871 11.793333 12.935484
## 1987  8.022581  9.810000 10.238710 13.150000 13.254839
## 1988  8.725806  9.883333 10.890323 12.253333 15.436667
## 1989  6.770968  8.486667  9.867742 12.876667 13.951613
## 1990  7.825806  9.166667 11.345161 12.656667 14.367742
plot(monthlyMinTempTS, col=c("blue"), ylab='Min temperature', xlab='Time', type='o',
main = "Time series plot of monthly minimum temperature in Melbourne.")

After converting the daily minimum temperature to a monthly mean max temperature, we can clearly see a strong seasonality pattern. Further Analysis will need to be performed.

# scatterPlot to Check correlations of changes in minimum daily temperature 

y = min_TempTS # read the daily min temp series into y 

x = zlag(min_TempTS) # generate first lag of the min temp series

head(y)
## Time Series:
## Start = c(1981, 1) 
## End = c(1981, 6) 
## Frequency = 365 
## [1] 20.7 17.9 18.8 14.6 15.8 15.8
head(x)
## [1]   NA 20.7 17.9 18.8 14.6 15.8
 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.7748702
plot(y[index],x[index],ylab='Min temperature series', xlab='The first lag of min temperature series', main = 'Scatter plot of min temperature series and its first lag')

As Scatter plot shows, the correlation value 0.774 proved that there is strong positive linear relationship between previous day’s min temperature on the next day’s min temperature

Investigating correlation of MonthlyMinTempTS

bb = monthlyMinTempTS
aa = zlag(monthlyMinTempTS)
index = 2: length(aa)
cor(bb[index], aa[index])
## [1] 0.8020674

Previous year temperature(Plotting the correlation)

We can see an upward trend and high correlation

plot(y = monthlyMinTempTS, x = zlag(monthlyMinTempTS), ylab = 'Min temp', xlab = 'Previous Year temp')

Trend analysis using regression approach

Time series plot of 1st difference of Minimum Monthly temperature levels

Although the general upward trend is resolved, we still have a seasonal effect pattern in time series plot.

plot(diff(monthlyMinTempTS),ylab='First Difference of Temperature',xlab='Time',main = "Figure-Time series plot of the first differences of Minimum Monthly temperature levels")

## ACF of first difference of Min Monthly temp levels

acf(as.vector(diff(monthlyMinTempTS)),lag.max=36,main="Figure-ACF of the first differences of Min monthly temp levels.")

PACF of first difference of Min Monthly temp levels

pacf(as.vector(diff(monthlyMinTempTS)),lag.max=36,main="Figure-PACF of the first differences of Min monthly temp levels.")

Although the general upward trend is resolved, we still have a seasonal effect pattern in time series plot. So, we consider taking a seasonal difference along with the ordinary difference

First and seasonal difference of Min Monthly temp levels

Time series plots of first and seasonal difference of Min Monthly temp levels

plot(diff(diff(monthlyMinTempTS),lag=12),xlab='Time',ylab='First and Seasonal Difference of monthly min temp levels',main= "Figure-Time series plot of the first and seasonal differences of monthly min temp levels.")

###ACF of the first and seasonal differences of Monthly min temp levels

acf(as.vector(diff(diff(monthlyMinTempTS),lag=12)),lag.max=36,ci.type='ma',main="Figure 12. ACF of the first and seasonal differences of")

###PACF of the first and seasonal differences of Monthly min temp levels

pacf(as.vector(diff(diff(monthlyMinTempTS),lag=12)),lag.max=36,ci.type='ma',main="Figure 12. PACF of the first and seasonal differences of monthly min temp levels")

After taking a seasonal difference, we nearly get rid of the seasonality in the series and very little autocorrelation remains in the series. This plot also suggests that a simple model which incorporates the lag 1 and lag 12 autocorrelations might be adequate. We will consider specifying an ARIMA(0,1,1)×(0,1,1)12 model to this series

##Residuals approach

###When D=1 ###Fit SARIMA(0,0,0)x(0,1,0)12 model and display time series plots of the residuals

m1.monthlyMinTempTS = arima(monthlyMinTempTS,order=c(0,0,0),seasonal=list(order=c(0,1,0), period=12))
res.m1 = residuals(m1.monthlyMinTempTS); 
plot(res.m1,xlab='Time',ylab='Residuals',main="m1-Time series plot of the residuals")

###ACF and PACF of residuals-m1

par(mfrow=c(1,2))
acf(res.m1, lag.max = 36, main = "m1-ACF of residuals")
pacf(res.m1, lag.max = 36, main = "m1-PACF of residuals")

When we focus on the autocorrelations at seasonal lags, there is no pattern implying the existence of a seasonal trend. However, the existence of an ordinary trend is obvious from the time series, ACF and PACF plots.

**(0,1,0)*(0,1,0)12 model**

To get rid of the ordinary trend seen in the residuals, we will fit the SARIMA(0,1,0)x(0,1,0)12 model

m2.monthlyMinTempTS = arima(monthlyMinTempTS,order=c(0,1,0),seasonal=list(order=c(0,1,0), period=12))
res.m2 = residuals(m2.monthlyMinTempTS); 
plot(res.m2,xlab='Time',ylab='Residuals',main="m2-Time series plot of the residuals")

###ACF and PACF of residuals-m2

par(mfrow=c(1,2))
acf(res.m2, lag.max = 36, main = "m2-ACF of residuals")
pacf(res.m2, lag.max = 36, main = "m2-PACF of residuals")

There is no evidence of an ordinary trend left in the latest residuals. The residuals of the model SARIMA(0,1,0)x(0,1,0)12 now includes SARMA and ARMA components.For the orders of SARMA component, we will consider the lags s,2s,3,…, i.e. 12, 24, 36, … These are shown as 1, 2, 3, … in both the ACF and PACF plots.

When we consider the lags 1, 2, 3, … of ACF plot, the correlation at lag 1 is significant, and in PACF there is a decreasing pattern at lags 1, 2, 3, … This implies the existence of an SMA(1) component.

**m3-(0,1,0)*(0,1,1)** Now, we will fit SARIMA(0,1,0)x(0,1,1)12 model and try to see if we get rid of the effect of the seasonal component on the residuals

m3.monthlyMinTempTS = arima(monthlyMinTempTS,order=c(0,1,0),seasonal=list(order=c(0,1,1), period=12))
res.m3 = residuals(m3.monthlyMinTempTS); 
plot(res.m3,xlab='Time',ylab='Residuals',main="m3-Time series plot of the residuals")

###ACF and PACF of residuals-m3

par(mfrow=c(1,2))
acf(res.m3, lag.max = 36, main = "m3-ACF of residuals")
pacf(res.m3, lag.max = 36, main = "m3-PACF of residuals")

At the seasonal lags, especially at the first seasonal lag in ACF, the autocorrelations became insignificant or slightly significant after adding the seasonal component. We will use the resulting ACF and PACF plots to specify the orders of ARMA component as there is no highly significant correlation at the lags s,2s,3s,….

The last ACF-PACF pair has significant pikes at the first lags implying an ARMA(1,1) model or due to the (insignificant) decrease of correlations in PACF, we would specify ARMA(0,1) for the residuals.

**m4-(0,1,1)*(0,1,1)

m4.monthlyMinTempTS = arima(monthlyMinTempTS,order=c(0,1,1),seasonal=list(order=c(0,1,1), period=12))
res.m4 = residuals(m4.monthlyMinTempTS); 
plot(res.m4,xlab='Time',ylab='Residuals',main="m4-Time series plot of the residuals")

###ACF and PACF of residuals-m4

par(mfrow=c(1,2))
acf(res.m4, lag.max = 36, main = "m4-ACF of residuals")
pacf(res.m4, lag.max = 36, main = "m4-PACF of residuals")

**m5-(1,1,1)*(0,1,1)

m5.monthlyMinTempTS = arima(monthlyMinTempTS,order=c(1,1,1),seasonal=list(order=c(0,1,1), period=12))
res.m5 = residuals(m5.monthlyMinTempTS); 
plot(res.m5,xlab='Time',ylab='Residuals',main="m5-Time series plot of the residuals")

###ACF and PACF of residuals-m5

par(mfrow=c(1,2))
acf(res.m5, lag.max = 36, main = "m5-ACF of residuals")
pacf(res.m5, lag.max = 36, main = "m5-PACF of residuals")

We got nearly white noise residual series for both of the SARIMA(0,1,1)x(0,1,1)12 and SARIMA(1,1,1)x(0,1,1)12 models. The residuals of SARIMA(0,1,1)x(0,1,1)12 model are closer to the white noise. So, we can conclude with the orders p=0, d=1, q=1, P=0, D=1, Q=1, and s=12 of SARIMA(p,d,q)x(P,D,Q)s model

EACF

we use EACF on the residuals of the non seasonary differencing model (i.e. res.m3), to see information about AR and MA parts lefts in residuals. The EACF table also clearly suggests a (0,1,1) model

eacf(res.m3)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o x x x o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x o x o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  o  o 
## 5 x x o x o o o o o o o  o  o  o 
## 6 o x x x o o o o o o o  o  o  o 
## 7 o x x o x o o o o o o  o  o  o

So, we can conclude from residual approach and EACF that the suitable models for this series are SARIMA(0,1,1)*(0,1,1)s=12 of SARIMA(p,d,q)x(P,D,Q)s model

Along with this, other candidate models are SARIMA(0,1,0) * (0,1,1)s=12 and (1,1,1)*(0,1,1)s=12

MODEL FITTING

PARAMETER ESTIMATION

  • Our final candidates for parameter estimation are – SARIMA(0,1,0)(0,1,1)S=12
    – SARIMA(1,1,1)    (0,1,1)S=12 – SARIMA(0,1,1)*(0,1,1)S=12

###**Fitting the SARIMA(0,1,0)*(0,1,1)S=12 model using ML method**

mf_010_ml = arima(monthlyMinTempTS,order=c(0,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "ML")
coeftest(mf_010_ml)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.785576   0.088072 -8.9197 < 2.2e-16 ***
## sma1 -0.801445   0.137944 -5.8099  6.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

**Fitting the same model above with lambda parameter to find AICc,BIC and AIC value

mf_010_lambda <- Arima(monthlyMinTempTS, order=c(0,1,0), seasonal=c(0,1,1),
  lambda=0)
mf_010_lambda
## Series: monthlyMinTempTS 
## ARIMA(0,1,0)(0,1,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          sma1
##       -0.9998
## s.e.   0.2219
## 
## sigma^2 estimated as 0.01344:  log likelihood=65.5
## AIC=-127   AICc=-126.88   BIC=-121.65

**Fitting the same SARIMA(0,1,0)*(0,1,1)S=12 model using CSS method**

mf_010_css = arima(monthlyMinTempTS,order=c(0,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "CSS")
coeftest(mf_010_css)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.787922   0.065622 -12.0070 < 2.2e-16 ***
## sma1 -0.586463   0.081213  -7.2213 5.151e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###**Fitting the SARIMA(1,1,1)*(0,1,1)S=12 model**

mf_111_ml = arima(monthlyMinTempTS,order=c(1,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "ML")
coeftest(mf_111_ml)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.079905   0.158437  0.5043     0.614    
## ma1  -0.836734   0.122602 -6.8248 8.804e-12 ***
## sma1 -0.824361   0.154570 -5.3333 9.646e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

**Fitting the same model above with lambda parameter to find AICc,BIC and AIC value

mf_111_lambda <- Arima(monthlyMinTempTS, order=c(1,1,1), seasonal=c(0,1,1),
  lambda=0)
mf_111_lambda
## Series: monthlyMinTempTS 
## ARIMA(1,1,1)(0,1,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ma1     sma1
##       -0.0083  -0.7410  -0.8670
## s.e.   0.2168   0.1982   0.1664
## 
## sigma^2 estimated as 0.009828:  log likelihood=88.59
## AIC=-169.17   AICc=-168.78   BIC=-158.48

**Fitting the same SARIMA(1,1,1)*(0,1,1)S=12 model using CSS method**

mf_111_css = arima(monthlyMinTempTS,order=c(1,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "CSS")
coeftest(mf_111_css)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.135634   0.140075 -0.9683    0.3329    
## ma1  -0.675582   0.117202 -5.7643 8.202e-09 ***
## sma1 -0.551512   0.084746 -6.5078 7.626e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###**Fitting the SARIMA(0,1,1)*(0,1,1)S=12 model**

mf_011_ml = arima(monthlyMinTempTS,order=c(0,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "ML")
coeftest(mf_011_ml)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.785576   0.088072 -8.9197 < 2.2e-16 ***
## sma1 -0.801445   0.137944 -5.8099  6.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

**Fitting the same model above with lambda parameter to find AICc,BIC and AIC value

mf_011_lambda <- Arima(monthlyMinTempTS, order=c(0,1,1), seasonal=c(0,1,1),
  lambda=0)
mf_011_lambda
## Series: monthlyMinTempTS 
## ARIMA(0,1,1)(0,1,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ma1     sma1
##       -0.7477  -0.8687
## s.e.   0.0894   0.1613
## 
## sigma^2 estimated as 0.009722:  log likelihood=88.59
## AIC=-171.17   AICc=-170.94   BIC=-163.15

**Fitting the same SARIMA(0,1,1)*(0,1,1)S=12 model using CSS method**

mf_011_css = arima(monthlyMinTempTS,order=c(0,1,1),seasonal=list(order=c(0,1,1),
                                      period=12),method = "CSS")
coeftest(mf_011_css)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.787922   0.065622 -12.0070 < 2.2e-16 ***
## sma1 -0.586463   0.081213  -7.2213 5.151e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###AIC values comparison The SARIMA model (0,1,1)*(0,1,1) has the lease AIC value and hence is considered to be the best fit

AIC(mf_010_lambda,mf_111_lambda,mf_011_lambda)
##               df       AIC
## mf_010_lambda  2 -126.9971
## mf_111_lambda  4 -169.1723
## mf_011_lambda  3 -171.1708

###BIC values comparison The SARIMA model (0,1,1)*(0,1,1) has the lease BIC value and hence is considered to be the best fit

BIC(mf_010_lambda,mf_111_lambda,mf_011_lambda)
##               df       BIC
## mf_010_lambda  2 -121.6515
## mf_111_lambda  4 -158.4809
## mf_011_lambda  3 -163.1523
  • **Because the standard error estimates for the SARIMA (0,1,1)*(0,1,1)s=12 model are relatively smaller than the corresponding parameter estimates, we think that the coefficient estimates are all highly significant**
  • **Because the SARIMA (0,1,1)*(0,1,1)s=12 model has the least AIC and BIC values, we declare this model as best fit model and to be used for forecasting**

DIAGNOSTIC CHECKING

plot(window(rstandard(mf_011_ml),start=c(1981,1)), ylab='Standardized Residuals',type='o', main="Residuals from the ARIMA(0,1,1)x(0,1,1)_12 model.")
abline(h=0)

Since we have some strange behaviour near 1983 , we will have to investigate this best model further for outliers

**plot the ACF of standardised residuals of our best model(0,1,1)*(0,1,1)s=12**

acf(as.vector(window(rstandard(mf_011_ml),start=c(1981,1))),lag.max=36,main="ACF of residuals from the ARIMA(0,1,1)x(0,1,1)_12 model.")

Except for the slightly significant autocorrelation at lag 28, there is not any sign of violation of the independence of residuals.

Function for residual analysis

residual.analysis <- function(model, std = TRUE){
  library(TSA)
  library(FitAR)
  if (std == TRUE){
    res.model = rstandard(model)
  }else{
    res.model = residuals(model)
  }
  par(mfrow=c(3,2))
  plot(res.model,type='o',ylab='Standardized residuals', main="Time series plot of standardized residuals")
  abline(h=0)
  hist(res.model,main="Histogram of standardized residuals")
  qqnorm(res.model,main="QQ plot of standardized residuals")
  qqline(res.model, col = 2)
  acf(res.model,main="ACF of standardized residuals")
  print(shapiro.test(res.model))
  k=0
  LBQPlot(res.model, lag.max = length(model$residuals)-1 , StartLag = k + 1, k = 0, SquaredQ = FALSE)
}

Please perform Box test,hist,qqnorm,qqline,shapiro-wilk test for all 3 models given below

**Residual Diagnostics - SARIMA(0,1,0)x(0,1,1)_12 Use this model name from above -> mf_010_ml

residual.analysis(model=mf_010_ml)
## 
##  Shapiro-Wilk normality test
## 
## data:  res.model
## W = 0.9922, p-value = 0.7394
## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple of
## shorter object length

SARIMA(0,1,0)x(0,1,1)_12 1. Time Series plots does not have any trend except for some intervention. 2. Histogram shows residuals normally distributed. 3. QQplot is normally distributed with 1 outlier. 4. ACF plot highlights the white noise presence. 5. Ljung Box test shows the presence of residuals.

**Residual Diagnostics - SARIMA(1,1,1)x(0,1,1)_12 ****Use this model name from above -> mf_111_ml**

residual.analysis(model=mf_111_ml)
## 
##  Shapiro-Wilk normality test
## 
## data:  res.model
## W = 0.99279, p-value = 0.7935
## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple of
## shorter object length

SARIMA(1,1,1)x(0,1,1)_12 1. Time Series plots does not have any trend except for some intervention. 2. Histogram is slightly skewed to the right. 3. QQplot is normally distributed except for some points not lying on the line. 4. ACF plot highlights the white noise presence. 5. Ljung Box test shows the presence of residuals.

**Residual Diagnostics - SARIMA(0,1,1)x(0,1,1)_12 Use this model name from above -> mf_011_ml

residual.analysis(model=mf_011_ml)
## 
##  Shapiro-Wilk normality test
## 
## data:  res.model
## W = 0.9922, p-value = 0.7394
## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple of
## shorter object length

SARIMA(0,1,1)x(0,1,1)_12 1. Time Series plots does not have any trend except for some intervention. 2. Histogram is slightly skewed to the right. 3. QQplot is normally distributed except for one point in outlier. 4. ACF plot highlights the white noise presence. 5. Ljung Box test shows the presence of residuals.

FORECAST

With SARIMA(0,1,1) * (0,1,1) being the best model with the S is equal to 12, below is the forcast for next 10 years.

sarima.for(monthlyMinTempTS,120,0,1,1,0,1,1,12)

## $pred
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1991 15.447837 15.650637 15.063298 12.787054 10.442440  7.810201  7.342003
## 1992 15.487655 15.690455 15.103116 12.826872 10.482258  7.850019  7.381821
## 1993 15.527473 15.730272 15.142934 12.866690 10.522075  7.889837  7.421638
## 1994 15.567290 15.770090 15.182752 12.906507 10.561893  7.929655  7.461456
## 1995 15.607108 15.809908 15.222570 12.946325 10.601711  7.969473  7.501274
## 1996 15.646926 15.849726 15.262388 12.986143 10.641529  8.009291  7.541092
## 1997 15.686744 15.889544 15.302206 13.025961 10.681347  8.049108  7.580910
## 1998 15.726562 15.929362 15.342023 13.065779 10.721165  8.088926  7.620728
## 1999 15.766380 15.969179 15.381841 13.105597 10.760982  8.128744  7.660545
## 2000 15.806198 16.008997 15.421659 13.145414 10.800800  8.168562  7.700363
##            Aug       Sep       Oct       Nov       Dec
## 1991  8.212052  9.443779 10.836433 12.950881 14.387385
## 1992  8.251870  9.483597 10.876251 12.990699 14.427203
## 1993  8.291688  9.523415 10.916069 13.030517 14.467020
## 1994  8.331506  9.563233 10.955887 13.070335 14.506838
## 1995  8.371324  9.603051 10.995705 13.110152 14.546656
## 1996  8.411141  9.642869 11.035523 13.149970 14.586474
## 1997  8.450959  9.682686 11.075341 13.189788 14.626292
## 1998  8.490777  9.722504 11.115158 13.229606 14.666110
## 1999  8.530595  9.762322 11.154976 13.269424 14.705927
## 2000  8.570413  9.802140 11.194794 13.309242 14.745745
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1991 1.024924 1.048196 1.070962 1.093255 1.115102 1.136529 1.157559 1.178215
## 1992 1.328409 1.354364 1.379831 1.404836 1.429405 1.453557 1.477315 1.500697
## 1993 1.659877 1.688194 1.716043 1.743447 1.770428 1.797003 1.823191 1.849008
## 1994 2.016801 2.047255 2.077263 2.106843 2.136014 2.164791 2.193191 2.221228
## 1995 2.397292 2.429717 2.461715 2.493302 2.524494 2.555306 2.585750 2.615840
## 1996 2.799858 2.834124 2.867980 2.901442 2.934522 2.967233 2.999587 3.031596
## 1997 3.223280 3.259281 3.294888 3.330115 3.364973 3.399474 3.433628 3.467445
## 1998 3.666531 3.704178 3.741447 3.778348 3.814893 3.851090 3.886950 3.922483
## 1999 4.128728 4.167947 4.206801 4.245299 4.283451 4.321266 4.358753 4.395921
## 2000 4.609103 4.649829 4.690201 4.730228 4.769920 4.809284 4.848329 4.887061
##           Sep      Oct      Nov      Dec
## 1991 1.198514 1.218475 1.238114 1.257447
## 1992 1.523721 1.546401 1.568754 1.590792
## 1993 1.874470 1.899590 1.924382 1.948859
## 1994 2.248916 2.276266 2.303292 2.330005
## 1995 2.645588 2.675005 2.704102 2.732889
## 1996 3.063271 3.094621 3.125657 3.156388
## 1997 3.500936 3.534110 3.566975 3.599540
## 1998 3.957697 3.992600 4.027200 4.061506
## 1999 4.432777 4.469329 4.505584 4.541550
## 2000 4.925489 4.963620 5.001460 5.039015

###Using auto.arima function To automatically fit the data we have and check the best model that could be fitted for our data using auto.arima function This result also comes as SARIMA(0,1,1)*(0,1,1)S=12

autoArimaModel = auto.arima(monthlyMinTempTS, d = 1)
autoArimaModel
## Series: monthlyMinTempTS 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.7856  -0.8014
## s.e.   0.0881   0.1379
## 
## sigma^2 estimated as 1.065:  log likelihood=-160.79
## AIC=327.58   AICc=327.81   BIC=335.6