Build.dat <- read.table("ApprovActiv.dat", header=T) ; attach(Build.dat)
App.ts <- ts(Approvals, start = c(1996,1), freq=4)
Act.ts <- ts(Activity, start = c(1996,1), freq=4)
ts.plot(App.ts, Act.ts, lty = c(1,3))
acf(ts.union(App.ts, Act.ts))
app.ran <- decompose(App.ts)$random
app.ran.ts <- window (app.ran, start = c(1996, 3) )
act.ran <- decompose (Act.ts)$random
act.ran.ts <- window (act.ran, start = c(1996, 3) )
T79 <- 1:10
Tdelt <- (1:100) / 10
Sales <- c(840,1470,2110,4000, 7590, 10950, 10530, 9470, 7790, 5890)
Cusales <- cumsum(Sales)
Bass.nls <- nls(Sales ~ M * ( ((P+Q)^2 / P) * exp(-(P+Q) * T79) ) /
(1+(Q/P)*exp(-(P+Q)*T79))^2, start = list(M=60630, P=0.03, Q=0.38))
summary(Bass.nls)
##
## Formula: Sales ~ M * (((P + Q)^2/P) * exp(-(P + Q) * T79))/(1 + (Q/P) *
## exp(-(P + Q) * T79))^2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## M 6.798e+04 3.128e+03 21.74 1.10e-07 ***
## P 6.594e-03 1.430e-03 4.61 0.00245 **
## Q 6.381e-01 4.140e-02 15.41 1.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 727.2 on 7 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 7.323e-06
Bcoef <- coef(Bass.nls)
m <- Bcoef[1]
p <- Bcoef[2]
q <- Bcoef[3]
ngete <- exp(-(p+q) * Tdelt)
Bpdf <- m * ( (p+q)^2 / p ) * ngete / (1 + (q/p) * ngete)^2
plot(Tdelt, Bpdf, xlab = "Year from 1979", ylab = "Sales per year", type='l')
points(T79, Sales)
Bcdf <- m * (1 - ngete)/(1 + (q/p)*ngete)
plot(Tdelt, Bcdf, xlab = "Year from 1979", ylab = "Cumulative sales", type='l')
points(T79, Cusales)
Motor.dat <- read.table("motororg.dat", header = T); attach(Motor.dat)
Comp.ts <- ts(complaints, start = c(1996, 1), freq = 12)
plot(Comp.ts, xlab = "Time / months", ylab = "Complaints")
wine.dat <- read.table("wine.dat", header = T) ; attach (wine.dat)
sweetw.ts <- ts(sweetw, start = c(1980,1), freq = 12)
plot(sweetw.ts, xlab= "Time (months)", ylab = "sales (1000 litres)")
sweetw.hw <- HoltWinters (sweetw.ts, seasonal = "mult")
sweetw.hw ; sweetw.hw$coef ; sweetw.hw$SSE
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = sweetw.ts, seasonal = "mult")
##
## Smoothing parameters:
## alpha: 0.4086698
## beta : 0
## gamma: 0.4929402
##
## Coefficients:
## [,1]
## a 285.6890314
## b 1.3509615
## s1 0.9498541
## s2 0.9767623
## s3 1.0275900
## s4 1.1991924
## s5 1.5463100
## s6 0.6730235
## s7 0.8925981
## s8 0.7557814
## s9 0.8227500
## s10 0.7241711
## s11 0.7434861
## s12 0.9472648
## a b s1 s2 s3 s4
## 285.6890314 1.3509615 0.9498541 0.9767623 1.0275900 1.1991924
## s5 s6 s7 s8 s9 s10
## 1.5463100 0.6730235 0.8925981 0.7557814 0.8227500 0.7241711
## s11 s12
## 0.7434861 0.9472648
## [1] 477693.9
sqrt(sweetw.hw$SSE/length(sweetw))
## [1] 50.54219
sd(sweetw)
## [1] 121.3908
plot (sweetw.hw$fitted)
plot (sweetw.hw)
HP.dat <- read.table("HP.txt", header = T) ; attach(HP.dat)
plot (as.ts(Price))
DP <- diff(Price) ; plot (as.ts(DP)) ; acf (DP)
mean(DP) + c(-2, 2) * sd(DP)/sqrt(length(DP))
## [1] 0.004378275 0.075353468
Global = scan("global.dat")
Global.ts = ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.ar <- ar(aggregate(Global.ts, FUN = mean), method = "mle")
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
mean(aggregate(Global.ts, FUN = mean))
## [1] -0.1382628
Global.ar$order
## [1] 4
Global.ar$ar
## [1] 0.58762026 0.01260254 0.11116731 0.26763656
acf(Global.ar$res[-(1:Global.ar$order)], lag = 50)
Global <- scan("global.dat")
Global.ts <- ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
temp <- window(Global.ts, start = 1970)
temp.lm <- lm(temp ~ time(temp))
coef(temp.lm)
## (Intercept) time(temp)
## -34.920409 0.017654
confint(temp.lm)
## 2.5 % 97.5 %
## (Intercept) -37.21001248 -32.63080554
## time(temp) 0.01650228 0.01880572
acf(resid(lm(temp ~ time(temp))))
library(nlme)
temp.gls <- gls(temp ~ time(temp), cor = corAR1(0.7))
confint(temp.gls)
## 2.5 % 97.5 %
## (Intercept) -39.80571681 -28.49658850
## time(temp) 0.01442274 0.02011148
Seas <- cycle(temp)
Time <- time(temp)
temp.lm <- lm(temp ~ 0 + Time + factor(Seas))
coef(temp.lm)
## Time factor(Seas)1 factor(Seas)2 factor(Seas)3 factor(Seas)4
## 0.01770758 -34.99726483 -34.98801824 -35.01002165 -35.01227506
## factor(Seas)5 factor(Seas)6 factor(Seas)7 factor(Seas)8 factor(Seas)9
## -35.03369513 -35.02505965 -35.02689640 -35.02476092 -35.03831988
## factor(Seas)10 factor(Seas)11 factor(Seas)12
## -35.05248996 -35.06557670 -35.04871900
new.t <- seq(2006, len = 2 * 12, by = 1/12)
alpha <- coef(temp.lm)[1]
beta <- rep(coef(temp.lm)[2:13], 2)
(alpha * new.t + beta)[1:4]
## factor(Seas)1 factor(Seas)2 factor(Seas)3 factor(Seas)4
## 0.5241458 0.5348681 0.5143403 0.5135625
new.dat <- data.frame(Time = new.t, Seas = rep(1:12, 2))
predict(temp.lm, new.dat)[1:24]
## 1 2 3 4 5 6 7
## 0.5241458 0.5348681 0.5143403 0.5135625 0.4936181 0.5037292 0.5033681
## 8 9 10 11 12 13 14
## 0.5069792 0.4948958 0.4822014 0.4705903 0.4889236 0.5418534 0.5525756
## 15 16 17 18 19 20 21
## 0.5320479 0.5312701 0.5113256 0.5214367 0.5210756 0.5246867 0.5126034
## 22 23 24
## 0.4999090 0.4882979 0.5066312
SIN <- COS <- matrix(nr = length(temp), nc = 6)
for (i in 1:6) {
COS[, i] <- cos(2 * pi * i * time(temp))
SIN[, i] <- sin(2 * pi * i * time(temp))
}
TIME <- (time(temp) - mean(time(temp)))/sd(time(temp))
mean(time(temp))
## [1] 1987.958
sd(time(temp))
## [1] 10.40433
temp.lm1 <- lm(temp ~ TIME + I(TIME^2) +
COS[,1] + SIN[,1] + COS[,2] + SIN[,2] +
COS[,3] + SIN[,3] + COS[,4] + SIN[,4] +
COS[,5] + SIN[,5] + COS[,6] + SIN[,6])
coef(temp.lm1)/sqrt(diag(vcov(temp.lm1)))
## (Intercept) TIME I(TIME^2) COS[, 1] SIN[, 1] COS[, 2]
## 18.2569117 30.2800742 1.2736899 0.7447906 2.3845512 1.2086419
## SIN[, 2] COS[, 3] SIN[, 3] COS[, 4] SIN[, 4] COS[, 5]
## 1.9310853 0.6448116 0.3971873 0.5468025 0.1681753 0.3169782
## SIN[, 5] COS[, 6] SIN[, 6]
## 0.3504607 -0.4498741 -0.6216650
temp.lm2 <- lm(temp ~ TIME + SIN[, 1] + SIN[, 2])
coef(temp.lm2)
## (Intercept) TIME SIN[, 1] SIN[, 2]
## 0.17501157 0.18409618 0.02041803 0.01615374
AIC(temp.lm)
## [1] -546.9577
AIC(temp.lm1)
## [1] -545.0491
AIC(temp.lm2)
## [1] -561.1779
plot(time(temp), resid(temp.lm2), type = "l")
abline(0, 0, col = "red")
acf(resid(temp.lm2))
pacf(resid(temp.lm2))
res.ar <- ar(resid(temp.lm2), method = "mle")
res.ar$ar
## [1] 0.4938189 0.3071598
sd(res.ar$res[-(1:2)])
## [1] 0.08373324
acf(res.ar$res[-(1:2)])
data(AirPassengers)
AP <- AirPassengers
plot(AP)
plot(log(AP))
SIN <- COS <- matrix(nr = length(AP), nc = 6)
for (i in 1:6) {
SIN[, i] <- sin(2 * pi * i * time(AP))
COS[, i] <- cos(2 * pi * i * time(AP))
}
TIME <- (time(AP) - mean(time(AP)))/sd(time(AP))
mean(time(AP))
## [1] 1954.958
sd(time(AP))
## [1] 3.476109
AP.lm1 <- lm(log(AP) ~ TIME + I(TIME^2) + I(TIME^3) + I(TIME^4) +
SIN[,1] + COS[,1] + SIN[,2] + COS[,2] + SIN[,3] + COS[,3] +
SIN[,4] + COS[,4] + SIN[,5] + COS[,5] + SIN[,6] + COS[,6])
coef(AP.lm1)/sqrt(diag(vcov(AP.lm1)))
## (Intercept) TIME I(TIME^2) I(TIME^3) I(TIME^4) SIN[, 1]
## 744.6853769 42.3818033 -4.1623354 -0.7514937 1.8725575 4.8678543
## COS[, 1] SIN[, 2] COS[, 2] SIN[, 3] COS[, 3] SIN[, 4]
## -26.0548660 10.3951486 10.0039542 -4.8436269 -1.5601366 -5.6662082
## COS[, 4] SIN[, 5] COS[, 5] SIN[, 6] COS[, 6]
## 1.9459153 -3.7662018 1.0259589 0.1500951 -0.5213860
AP.lm2 <- lm(log(AP) ~ TIME + I(TIME^2) + SIN[,1] + COS[,1] +
SIN[,2] + COS[,2] + SIN[,3] + SIN[,4] + COS[,4] + SIN[,5])
coef(AP.lm2)/sqrt(diag(vcov(AP.lm2)))
## (Intercept) TIME I(TIME^2) SIN[, 1] COS[, 1] SIN[, 2]
## 922.625445 103.520628 -8.235569 4.916041 -25.813114 10.355319
## COS[, 2] SIN[, 3] SIN[, 4] COS[, 4] SIN[, 5]
## 9.961335 -4.789997 -5.612235 1.949391 -3.730666
AIC(AP.lm1)
## [1] -447.9475
AIC(AP.lm2)
## [1] -451.0749
acf(resid(AP.lm2))
AP.gls <- gls(log(AP) ~ TIME + I(TIME^2) + SIN[,1] + COS[,1] +
SIN[,2] + COS[,2] + SIN[,3] + SIN[,4] + COS[,4] + SIN[,5],
cor = corAR1(0.6))
coef(AP.gls)/sqrt(diag(vcov(AP.gls)))
## (Intercept) TIME I(TIME^2) SIN[, 1] COS[, 1] SIN[, 2]
## 398.841577 45.850520 -3.651287 3.299713 -18.177844 11.768934
## COS[, 2] SIN[, 3] SIN[, 4] COS[, 4] SIN[, 5]
## 11.432678 -7.630660 -10.749938 3.571395 -7.920777
AP.ar <- ar(resid(AP.lm2), order = 1, method = "mle")
AP.ar$ar
## [1] 0.6410685
acf(AP.ar$res[-1])
new.t <- time(ts(start = 1961, end = c(1970, 12), fr = 12))
TIME <- (new.t - mean(time(AP)))/sd(time(AP))
SIN <- COS <- matrix(nr = length(new.t), nc = 6)
for (i in 1:6) {
COS[, i] <- cos(2 * pi * i * new.t)
SIN[, i] <- sin(2 * pi * i * new.t)
}
SIN <- SIN[, -6]
new.dat <- data.frame(TIME = as.vector(TIME), SIN = SIN,
COS = COS)
AP.pred.ts <- exp(ts(predict(AP.lm2, new.dat), st = 1961, fr = 12))
ts.plot(log(AP), log(AP.pred.ts), lty = 1:2)
ts.plot(AP, AP.pred.ts, lty = 1:2)
summary(AP.lm2)$r.sq
## [1] 0.9888318
sigma <- summary(AP.lm2)$sigma
lognorm.correction.factor <- exp((1/2) * sigma^2)
empirical.correction.factor <- mean(exp(resid(AP.lm2)))
lognorm.correction.factor
## [1] 1.001171
empirical.correction.factor
## [1] 1.00108
AP.pred.ts <- AP.pred.ts * empirical.correction.factor
x <- read.table("pounds_nz.dat", header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
##
## Call:
## arima(x = x.ts, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 1.000 2.8329
## s.e. 0.072 0.0646
##
## sigma^2 estimated as 0.04172: log likelihood = 4.76, aic = -3.53
acf(x.ma$res[-1])
set.seed(1)
x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1)))
## ar1 ma1 intercept
## -0.596966371 0.502703368 -0.006571345
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.arma <- arima(x.ts, order = c(1, 0, 1))
AIC(x.ma)
## [1] -3.526895
AIC(x.ar)
## [1] -37.40417
AIC(x.arma)
## [1] -42.27357
x.arma
##
## Call:
## arima(x = x.ts, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.8925 0.5319 2.9597
## s.e. 0.0759 0.2021 0.2435
##
## sigma^2 estimated as 0.01505: log likelihood = 25.14, aic = -42.27
acf(resid(x.arma))
CBE <- read.table("cbe.dat", header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
Time <- 1:length(Elec.ts)
Imth <- cycle(Elec.ts)
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
acf(resid(Elec.lm))
best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:2) for (j in 0:2) {
fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0,
j)))
if (fit.aic < best.aic) {
best.order <- c(i, 0, j)
best.arma <- arima(resid(Elec.lm), order = best.order)
best.aic <- fit.aic
}
}
best.order
## [1] 2 0 0
acf(resid(best.arma))
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12, 3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991, freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)
wave.dat <- read.table("wave.dat", header = T)
attach (wave.dat)
layout(1:3)
plot (as.ts(waveht), ylab = 'Wave height')
acf (waveht)
pacf (waveht)
wave.arma <- arima(waveht, order = c(4,0,4))
acf (wave.arma$res[-(1:4)])
pacf (wave.arma$res[-(1:4)])
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')
CBE <- read.table("cbe.dat", he = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
layout(c(1, 1, 2, 3))
plot(Elec.ts)
plot(diff(Elec.ts))
plot(diff(log(Elec.ts)))
set.seed(1)
x <- w <- rnorm(1000)
for (i in 3:1000) x[i] <- 0.5 * x[i - 1] + x[i - 1] - 0.5 *
x[i - 2] + w[i] + 0.3 * w[i - 1]
arima(x, order = c(1, 1, 1))
##
## Call:
## arima(x = x, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.4235 0.3308
## s.e. 0.0433 0.0450
##
## sigma^2 estimated as 1.067: log likelihood = -1450.13, aic = 2906.26
x <- arima.sim(model = list(order = c(1, 1, 1), ar = 0.5, ma = 0.3), n = 1000)
arima(x, order = c(1, 1, 1))
##
## Call:
## arima(x = x, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.5567 0.2502
## s.e. 0.0372 0.0437
##
## sigma^2 estimated as 1.079: log likelihood = -1457.45, aic = 2920.91
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Beer.ima <- arima(Beer.ts, order = c(0, 1, 1))
Beer.ima
##
## Call:
## arima(x = Beer.ts, order = c(0, 1, 1))
##
## Coefficients:
## ma1
## -0.3334
## s.e. 0.0558
##
## sigma^2 estimated as 360.4: log likelihood = -1723.27, aic = 3450.53
Beer.1991 <- predict(Beer.ima, n.ahead = 12)
sum(Beer.1991$pred)
## [1] 2365.412
AIC (arima(log(Elec.ts), order = c(1,1,0), seas = list(order = c(1,0,0), 12)))
## [1] -1764.741
AIC (arima(log(Elec.ts), order = c(0,1,1), seas = list(order = c(0,0,1), 12)))
## [1] -1361.586
library(MASS)
data(SP500)
plot(SP500, type = 'l')
acf(SP500)
acf((SP500 - mean(SP500))^2)
set.seed(1)
alpha0 <- 0.1
alpha1 <- 0.4
beta1 <- 0.2
w <- rnorm(10000)
a <- rep(0, 10000)
h <- rep(0, 10000)
for (i in 2:10000) {
h[i] <- alpha0 + alpha1 * (a[i - 1]^2) + beta1 * h[i - 1]
a[i] <- w[i] * sqrt(h[i])
}
acf(a)
acf(a^2)
library(tseries)
a.garch <- garch(a, grad = "numerical", trace = FALSE)
confint(a.garch)
## 2.5 % 97.5 %
## a0 0.0882393 0.1092903
## a1 0.3307897 0.4023932
## b1 0.1928344 0.2954660
sp.garch <- garch(SP500, trace = F)
sp.res <- sp.garch$res[-1]
acf(sp.res)
acf(sp.res^2)
get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1))
{
best.aic <- 1e8
n <- length(x.ts)
for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
{
fit <- arima(x.ts, order = c(p,d,q),
seas = list(order = c(P,D,Q),
frequency(x.ts)), method = "CSS")
fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
if (fit.aic < best.aic)
{
best.aic <- fit.aic
best.fit <- fit
best.model <- c(p,d,q,P,D,Q)
}
}
list(best.aic, best.fit, best.model)
}
stemp <- scan("stemp.dat")
stemp.ts <- ts(stemp, start = 1850, freq = 12)
plot(stemp.ts)
stemp.best <- get.best.arima(stemp.ts, maxord = rep(2,6))
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D,
## Q), : possible convergence problem: optim gave code = 1
stemp.best[[3]]
## [1] 1 1 2 2 0 1
stemp.arima <- arima(stemp.ts, order = c(1,1,2),
seas = list(order = c(2,0,1), 12))
t(confint(stemp.arima) )
## ar1 ma1 ma2 sar1 sar2 sma1
## 2.5 % 0.8317386 -1.447401 0.3256697 0.8576773 -0.02501929 -0.9690566
## 97.5 % 0.9127950 -1.312552 0.4530477 1.0041425 0.07413469 -0.8507002
stemp.res <- resid(stemp.arima)
layout(1:2)
acf(stemp.res)
acf(stemp.res^2)
stemp.garch <- garch(stemp.res, trace = F)
t(confint(stemp.garch))
## a0 a1 b1
## 2.5 % 9.710626e-06 0.03264244 0.9257899
## 97.5 % 1.462336e-04 0.06455859 0.9635381
stemp.garch.res <- resid(stemp.garch)[-1]
acf(stemp.garch.res)
acf(stemp.garch.res^2)
xrates <- read.table("us_rates.dat", header = T)
xrates[1:3, ]
## UK NZ EU
## 1 0.55834 1.5168 0.79419
## 2 0.55288 1.4884 0.78860
## 3 0.54803 1.4885 0.78284
acf( diff(xrates$UK) )
acf( diff(xrates$EU) )
plot(xrates$UK, xrates$EU, pch = 4)
cor(xrates$UK, xrates$EU)
## [1] 0.9462121
library(tseries)
adf.test(x)
##
## Augmented Dickey-Fuller Test
##
## data: x
## Dickey-Fuller = -2.7401, Lag order = 9, p-value = 0.265
## alternative hypothesis: stationary
pp.test(xrates$UK)
##
## Phillips-Perron Unit Root Test
##
## data: xrates$UK
## Dickey-Fuller Z(alpha) = -10.554, Truncation lag parameter = 7,
## p-value = 0.521
## alternative hypothesis: stationary
pp.test(xrates$EU)
##
## Phillips-Perron Unit Root Test
##
## data: xrates$EU
## Dickey-Fuller Z(alpha) = -6.812, Truncation lag parameter = 7,
## p-value = 0.7297
## alternative hypothesis: stationary
x <- y <- mu <- rep(0, 1000)
for (i in 2:1000) mu[i] <- mu[i - 1] + rnorm(1)
x <- mu + rnorm(1000)
y <- mu + rnorm(1000)
adf.test(x)$p.value
## [1] 0.02580509
adf.test(y)$p.value
## [1] 0.03060587
po.test(cbind(x, y))
## Warning in po.test(cbind(x, y)): p-value smaller than printed p-value
##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(x, y)
## Phillips-Ouliaris demeaned = -1004.9, Truncation lag parameter =
## 9, p-value = 0.01
po.test(cbind(xrates$UK, xrates$EU))
##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(xrates$UK, xrates$EU)
## Phillips-Ouliaris demeaned = -21.662, Truncation lag parameter =
## 10, p-value = 0.04118
ukeu.lm <- lm(xrates$UK ~ xrates$EU)
ukeu.res <- resid(ukeu.lm)
ukeu.res.ar <- ar(ukeu.res)
ukeu.res.ar$order
## [1] 3
AIC(arima(ukeu.res, order = c(3, 0, 0)))
## [1] -9886.26
AIC(arima(ukeu.res, order = c(2, 0, 0)))
## [1] -9886.157
AIC(arima(ukeu.res, order = c(1, 0, 0)))
## [1] -9879.82
AIC(arima(ukeu.res, order = c(1, 1, 0)))
## [1] -9875.723
library(tseries)
data(USeconomic)
US.ar <- ar(cbind(GNP, M1), method="ols", dmean=T, intercept=F)
US.ar$ar
## , , GNP
##
## GNP M1
## 1 1.271812104 -0.03383385
## 2 -0.004229937 0.06353801
## 3 -0.267154022 -0.02858942
##
## , , M1
##
## GNP M1
## 1 1.1674655 1.5876695
## 2 -0.6941813 -0.4838919
## 3 -0.5103451 -0.1294549
acf(US.ar$res[-c(1:3), 1])
acf(US.ar$res[-c(1:3), 2])
library(vars)
## Loading required package: strucchange
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
US.var <- VAR(cbind(GNP, M1), p = 3, type = "trend")
coef(US.var)
## $GNP
## Estimate Std. Error t value Pr(>|t|)
## GNP.l1 1.075374027 0.08843033 12.1606917 5.478655e-23
## M1.l1 1.036148232 0.41028680 2.5254243 1.279623e-02
## GNP.l2 -0.006780841 0.13281676 -0.0510541 9.593633e-01
## M1.l2 -0.300378560 0.75433064 -0.3982054 6.911527e-01
## GNP.l3 -0.127236772 0.08508606 -1.4953891 1.373134e-01
## M1.l3 -0.563697120 0.44568891 -1.2647771 2.082856e-01
## trend 1.035029562 0.43519492 2.3783126 1.889452e-02
##
## $M1
## Estimate Std. Error t value Pr(>|t|)
## GNP.l1 -0.04391500 0.01910790 -2.2982638 2.319447e-02
## M1.l1 1.59230681 0.08865419 17.9608749 1.505956e-36
## GNP.l2 0.06164770 0.02869886 2.1480889 3.362031e-02
## M1.l2 -0.48913808 0.16299470 -3.0009448 3.245326e-03
## GNP.l3 -0.01754739 0.01838528 -0.9544264 3.416953e-01
## M1.l3 -0.10410852 0.09630383 -1.0810424 2.817430e-01
## trend 0.01159781 0.09403630 0.1233333 9.020397e-01
US.var <- VAR(cbind(GNP, M1), p = 2, type = "trend")
coef(US.var)
## $GNP
## Estimate Std. Error t value Pr(>|t|)
## GNP.l1 1.1413956 0.08448544 13.509967 1.831652e-26
## M1.l1 1.3304479 0.33914181 3.922984 1.414205e-04
## GNP.l2 -0.1996175 0.08230909 -2.425218 1.668345e-02
## M1.l2 -1.1566875 0.34879323 -3.316255 1.185205e-03
## trend 1.0319351 0.42297164 2.439726 1.605720e-02
##
## $M1
## Estimate Std. Error t value Pr(>|t|)
## GNP.l1 -0.03371829 0.01810619 -1.86225184 6.484268e-02
## M1.l1 1.64898451 0.07268195 22.68767536 7.330886e-47
## GNP.l2 0.03419264 0.01763978 1.93838296 5.476075e-02
## M1.l2 -0.65016206 0.07475036 -8.69777827 1.345903e-14
## trend 0.00654324 0.09064764 0.07218324 9.425679e-01
acf(resid(US.var)[, 1])
acf(resid(US.var)[, 2])
US.pred <- predict(US.var, n.ahead = 4)
US.pred
## $GNP
## fcst lower upper CI
## [1,] 3957.507 3911.310 4003.704 46.19721
## [2,] 3986.346 3913.707 4058.986 72.63959
## [3,] 4014.460 3921.493 4107.426 92.96630
## [4,] 4042.940 3933.009 4152.872 109.93141
##
## $M1
## fcst lower upper CI
## [1,] 630.8020 620.9014 640.7026 9.900588
## [2,] 631.8559 612.8800 650.8318 18.975878
## [3,] 633.6114 606.1469 661.0759 27.464512
## [4,] 635.8656 600.7831 670.9482 35.082583
GNP.pred <- ts(US.pred$fcst$GNP[, 1], st = 1988, fr = 4)
M1.pred <- ts(US.pred$fcst$M1[, 1], st = 1988, fr = 4)
ts.plot(cbind(window(GNP, start = 1981), GNP.pred), lty = 1:2)
ts.plot(cbind(window(M1, start = 1981), M1.pred), lty = 1:2)
