Practical Time Series Analysis - Week 4 Labs
Ken Wood
12/20/2021
Partial Auto-Correlation Function (PACF)
phi.1 = .9; phi.2 = -.6;
par(mfrow=c(3,1))
data.ts = arima.sim(n = 500, list(ar = c(phi.1, phi.2)))
plot(data.ts,
main =
paste("Autoregressive Process with phi1=", phi.1," phi2=",phi.2 ) )
acf(data.ts) #ACF
acf(data.ts, type="partial") #partial ACF
Beveridge Wheat Price Data Analysis
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoodata(bev)
plot(bev,ylab='Price',main='Beveridge Wheat Price Data')
bev.MA <- filter(bev,rep(1/31,31),sides = 2)
lines(bev.MA,col='red')
par(mfrow=c(3,1))
y <- bev/bev.MA
plot(y,ylab='Scaled Price', main='Transformed Beveridge Wheat Price Data')
acf(na.omit(y),main='ACF of Transformed Beveridge Wheat Price Data')
acf(na.omit(y),type='partial',main='PACF of Transformed Beveridge Wheat Price Data')
Body Fat Data
library(isdals)
data(bodyfat)
attach(bodyfat) # `bodyfat` columns can now be directly accessed.
pairs(bodyfat)
cor(bodyfat)## Fat Triceps Thigh Midarm
## Fat 1.0000000 0.8432654 0.8780896 0.1424440
## Triceps 0.8432654 1.0000000 0.9238425 0.4577772
## Thigh 0.8780896 0.9238425 1.0000000 0.0846675
## Midarm 0.1424440 0.4577772 0.0846675 1.0000000# Create a linear model to predict `Fat`
Fat.hat <- predict(lm(Fat~Thigh))
Triceps.hat <- predict(lm(Triceps~Thigh))
cor(Fat-Fat.hat,Triceps-Triceps.hat) # correlation of residuals## [1] 0.1749822library(ppcor)## Loading required package: MASSpcor(cbind(Fat,Triceps,Thigh)) # we get the same answer using the partial correlation function## $estimate
## Fat Triceps Thigh
## Fat 1.0000000 0.1749822 0.4814109
## Triceps 0.1749822 1.0000000 0.7130120
## Thigh 0.4814109 0.7130120 1.0000000
##
## $p.value
## Fat Triceps Thigh
## Fat 0.00000000 0.4736789763 0.0368987227
## Triceps 0.47367898 0.0000000000 0.0006109801
## Thigh 0.03689872 0.0006109801 0.0000000000
##
## $statistic
## Fat Triceps Thigh
## Fat 0.0000000 0.7327755 2.264597
## Triceps 0.7327755 0.0000000 4.192849
## Thigh 2.2645969 4.1928494 0.000000
##
## $n
## [1] 20
##
## $gp
## [1] 1
##
## $method
## [1] "pearson"Fat.hat <- predict(lm(Fat~Thigh+Midarm))
Triceps.hat <- predict(lm(Triceps~Thigh+Midarm))
cor(Fat-Fat.hat,Triceps-Triceps.hat) # correlation of residuals## [1] 0.33815pcor(bodyfat) # we get the same answer using the partial correlation function## $estimate
## Fat Triceps Thigh Midarm
## Fat 1.0000000 0.3381500 -0.2665991 -0.3240520
## Triceps 0.3381500 1.0000000 0.9963725 0.9955918
## Thigh -0.2665991 0.9963725 1.0000000 -0.9926612
## Midarm -0.3240520 0.9955918 -0.9926612 1.0000000
##
## $p.value
## Fat Triceps Thigh Midarm
## Fat 0.0000000 1.699111e-01 2.848944e-01 1.895628e-01
## Triceps 0.1699111 0.000000e+00 1.490492e-18 7.071386e-18
## Thigh 0.2848944 1.490492e-18 0.000000e+00 4.134178e-16
## Midarm 0.1895628 7.071386e-18 4.134178e-16 0.000000e+00
##
## $statistic
## Fat Triceps Thigh Midarm
## Fat 0.000000 1.437266 -1.106441 -1.370142
## Triceps 1.437266 0.000000 46.833317 42.459264
## Thigh -1.106441 46.833317 0.000000 -32.834695
## Midarm -1.370142 42.459264 -32.834695 0.000000
##
## $n
## [1] 20
##
## $gp
## [1] 2
##
## $method
## [1] "pearson"Lake Huron data
par( mfrow=c(3,1) )
plot(LakeHuron)
acf(LakeHuron)
acf(LakeHuron, type="partial")
Yule Walker Estimation Simulation of AR(2) process
\[\begin{align} x_t=\phi_1x_{t-1}+\phi_2x_{t-2}+z_t \newline z_t \sim N(0, \sigma^2) \end{align}\]
# set seed a common number, so we can reproduce the same datasets
set.seed(2017)# model parameters (we will estimate them)
sigma=4
phi=NULL
phi[1:2]=c(1/3,1/2)
phi## [1] 0.3333333 0.5000000# number of data points
n=10000# Simulate AR process
ar.process=arima.sim(n,model=list(ar=c(1/3,1/2)), sd=4)
ar.process[1:5]## [1] 4.087685 5.598492 3.019295 2.442354 5.398302# find and name 2nd and 3rd sample autocorrelation
r=NULL
r[1:2]=acf(ar.process, plot=F)$acf[2:3]
r## [1] 0.6814103 0.7255825# matrix R
R=matrix(1,2,2) # matrix of dimension 2 by 2, with entries all 1's.
R## [,1] [,2]
## [1,] 1 1
## [2,] 1 1# edit R
R[1,2]=r[1] # only diagonal entries are edited
R[2,1]=r[1] # only diagonal entries are edited
R## [,1] [,2]
## [1,] 1.0000000 0.6814103
## [2,] 0.6814103 1.0000000# b-column vector on the right (transformed from `r` above)
b=matrix(r,nrow=2,ncol=1)# b- column vector with no entries
b## [,1]
## [1,] 0.6814103
## [2,] 0.7255825\(\hat{b}=\hat{R}\hat{\phi}\). Solving for \(\hat{\phi}\), we obtain the \(\hat{\phi}=\hat{R}^{-1}\hat{b}\) vector
phi.hat=matrix(c(solve(R,b)[1,1], solve(R,b)[2,1]),2,1)
phi.hat## [,1]
## [1,] 0.3490720
## [2,] 0.4877212# variance estimation
c0=acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 16.37169# plot time series, along with acf, pacf
par(mfrow=c(3,1))
plot(ar.process, main='Simulated AR(2)')
acf(ar.process, main='ACF')
pacf(ar.process, main='PACF')
Yule Walker Estimation of model patrameters of an AR(3) simulation
\[\begin{align} x_t=\phi_1x_{t-1}+\phi_2x_{t-2}+\phi_3x_{t-3}+z_t \newline z_t \sim N(0, \sigma^2) \end{align}\]
# set seed a common number, so we can reproduce the same datasets
set.seed(2017)# model parameters (we will estimate them)
sigma=4
phi=NULL
phi[1:3]=c(1/3,1/2,7/100)
phi## [1] 0.3333333 0.5000000 0.0700000# number of data points
n=10000# Simulate AR(3) process
ar3.process=arima.sim(n,model=list(ar=c(1/3,1/2, 7/100)), sd=4)# find and name 2nd, 3rd and 4th sample autocorrelation
r=NULL
r[1:3]=acf(ar3.process, plot=F)$acf[2:4]
r## [1] 0.7894308 0.8167192 0.7367485# matrix R
R=matrix(1,3,3) # matrix of dimension 3 by 3, with entries all 1's.
R## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 1 1
## [3,] 1 1 1# Populate the R matrix
R[1,2]=r[1]
R[1,3]=r[2]
R[2,1]=r[1]
R[2,3]=r[1]
R[3,1]=r[2]
R[3,2]=r[1]
R## [,1] [,2] [,3]
## [1,] 1.0000000 0.7894308 0.8167192
## [2,] 0.7894308 1.0000000 0.7894308
## [3,] 0.8167192 0.7894308 1.0000000# b-column vector on the right
b=matrix(,3,1) # b- column vector with no entries
b[1,1]=r[1]
b[2,1]=r[2]
b[3,1]=r[3]
b## [,1]
## [1,] 0.7894308
## [2,] 0.8167192
## [3,] 0.7367485# solve Rx=b and find phi's
phi.hat=solve(R,b)
phi.hat## [,1]
## [1,] 0.3512298
## [2,] 0.4890902
## [3,] 0.0637896# sigme estimation
c0=acf(ar3.process, type='covariance', plot=F)$acf[1]
var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 16.38984Plots
par(mfrow=c(3,1))
plot(ar3.process, main='Simulated AR(3)')
acf(ar3.process, main='ACF')
pacf(ar3.process, main='PACF')
Modeling recruitment time series from ‘astsa’ package as an AR process
library(astsa)
my.data=rec
# Plot rec
plot(rec, main='Recruitment time series', col='blue', lwd=3)
# subtract mean to get a time series with mean zero
ar.process=my.data-mean(my.data)# ACF and PACF of the process
par(mfrow=c(2,1))
acf(ar.process, main='Recruitment', col='red', lwd=3)
pacf(ar.process, main='Recruitment', col='green', lwd=3)
# order
p=2
# sample autocorreleation function r
r=NULL
r[1:p]=acf(ar.process, plot=F)$acf[2:(p+1)]
cat('r=',r,'\n')## r= 0.9218042 0.7829182# matrix R
R=matrix(1,p,p) # matrix of dimension 2 by 2, with entries all 1's.
# define non-diagonal entires of R
for(i in 1:p){
for(j in 1:p){
if(i!=j)
R[i,j]=r[abs(i-j)]
}
}
R## [,1] [,2]
## [1,] 1.0000000 0.9218042
## [2,] 0.9218042 1.0000000# b-column vector on the right
b=NULL
b=matrix(r,p,1)# b- column vector with no entries
b## [,1]
## [1,] 0.9218042
## [2,] 0.7829182# solve(R,b) solves Rx=b, and gives x=R^(-1)b vector
phi.hat=NULL
phi.hat=solve(R,b)[,1]
phi.hat## [1] 1.3315874 -0.4445447#variance estimation using Yule-Walker Estimator
c0=acf(ar.process, type='covariance', plot=F)$acf[1]
c0## [1] 780.991var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 94.17131# constant term in the model
phi0.hat=mean(my.data)*(1-sum(phi.hat))
phi0.hat## [1] 7.033036cat("Constant:", phi0.hat," Coefficients:", phi.hat, " and Variance:", var.hat, '\n')## Constant: 7.033036 Coefficients: 1.331587 -0.4445447 and Variance: 94.17131Johnson & Johnson quarterly earnings per share data
# Time plot for Johnson&Johnson
plot(JohnsonJohnson, main='Johnson&Johnosn earnings per share', col='blue', lwd=3)
# log-return of Johnson&Johnson
jj.log.return=diff(log(JohnsonJohnson))
jj.log.return.mean.zero=jj.log.return-mean(jj.log.return)
# Plots for log-returns
par(mfrow=c(3,1))
plot(jj.log.return.mean.zero, main='Log-return (mean zero) of Johnson&Johnosn earnings per share')
acf(jj.log.return.mean.zero, main='ACF')
pacf(jj.log.return.mean.zero, main='PACF')
# Order
p=4
# sample autocorreleation function r
r=NULL
r[1:p]=acf(jj.log.return.mean.zero, plot=F)$acf[2:(p+1)]
r## [1] -0.50681760 0.06710084 -0.40283604 0.73144780# matrix R
R=matrix(1,p,p) # matrix of dimension 4 by 4, with entries all 1's.
# define non-diagonal entires of R
for(i in 1:p){
for(j in 1:p){
if(i!=j)
R[i,j]=r[abs(i-j)]
}
}
R## [,1] [,2] [,3] [,4]
## [1,] 1.00000000 -0.50681760 0.06710084 -0.40283604
## [2,] -0.50681760 1.00000000 -0.50681760 0.06710084
## [3,] 0.06710084 -0.50681760 1.00000000 -0.50681760
## [4,] -0.40283604 0.06710084 -0.50681760 1.00000000# b-column vector on the right
b=matrix(r,p,1)# b- column vector with no entries
b## [,1]
## [1,] -0.50681760
## [2,] 0.06710084
## [3,] -0.40283604
## [4,] 0.73144780phi.hat=solve(R,b)[,1]
phi.hat## [1] -0.6293492 -0.5171526 -0.4883374 0.2651266# Variance estimation using Yule-Walker Estimator
c0=acf(jj.log.return.mean.zero, type='covariance', plot=F)$acf[1]
c0## [1] 0.04365692var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 0.01419242# Constant term in the model
phi0.hat=mean(jj.log.return)*(1-sum(phi.hat))
phi0.hat## [1] 0.079781cat("Constant:", phi0.hat," Coeffcinets:", phi.hat, " and Variance:", var.hat, '\n')## Constant: 0.079781 Coeffcinets: -0.6293492 -0.5171526 -0.4883374 0.2651266 and Variance: 0.01419242