varve = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/varve.dat"))
varve = ts(scan("varve.dat"))

## 'scan' for univariate series, 'read.table' for multivariate

## Source the supplied functions
source("http://www.stat.pitt.edu/stoffer/tsa2/Rcode/itall.R")

## The first differences of the logs look fairly stationary
y = diff(log(varve))
par(mfrow=c(3,1))
plot(varve)
plot(log(varve))
plot(y)


## Fit an MA(1):
fit1.y = sarima(y, p = 0, d = 0, q = 1)
fit1.y

resids = fit1.y$fit$resid
pacf(resids)

## Fit an ARMA(1,1):
fit2.y = sarima(y, p = 1, d = 0, q = 1)
fit2.y

## Use the better model to forecast:

forecasts.varve = sarima.for(log(varve), nahead = 100, p=1, d=1, q=1)
forecasts.varve