Practice Problems 4

library(Matrix)
A <- read.table("/Users/YaseenAbdulridha/Downloads/3_06.csv")
x <- A$x
y <- A$y
z <- A$z


log_mvnpdf = function( v, mu, A){
  n<- length(v)
  u <- v-mu
  ss <- sum(solve(A,u)*u)
  det <- 2*sum( log( diag( chol( A )))) 
  logL <- .5*( ss + det + n*log(2*pi) )
}

AR(2)

xlogL = function( theta ) {
  n     <- length(x)
  mu    <- theta[1]
  phi1  <- theta[2]
  phi2  <- theta[3]
  sigma <- theta[4]
  ii <- 3:n
  logL <- -sum( dnorm( x[ii], mu + phi1*x[ii-1] + phi2*x[ii-2], sigma, log=TRUE) )
}
optim( c(0,0,0,1), xlogL, method="BFGS" )

## $par
## [1] -0.01628134  0.95735919 -0.22430406  0.97910072
## 
## $value
## [1] 1395.021
## 
## $counts
## function gradient 
##       62       17 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

arima( x, order=c(2,0,0) )

## 
## Call:
## arima(x = x, order = c(2, 0, 0))
## 
## Coefficients:
##          ar1      ar2  intercept
##       0.9571  -0.2241    -0.0690
## s.e.  0.0308   0.0309     0.1157
## 
## sigma^2 estimated as 0.9578:  log likelihood = -1397.91,  aic = 2803.82

MA(2)

ylogL = function( theta) {
  n     <- length(y)
  mu    <- theta[1]
  1
  rho   <- ARMAacf( ma=theta[2:3] )
  gamma <- rho * theta[4]^2 * ( 1 + sum(theta[2:3]^2 ))
  A <- bandSparse( n,
                  n,
                  k=0:2,
                  diagonals=cbind( rep(gamma[1],n), rep(gamma[2],n), rep(gamma[3],n) ),
                  symmetric=TRUE)
  logL <- log_mvnpdf( y, mu, A )
}
optim( c(0,0,0,1), ylogL, method="BFGS" )

## $par
## [1] -0.02507634  0.71874704 -0.16294877  1.09991286
## 
## $value
## [1] 1514.871
## 
## $counts
## function gradient 
##       33       10 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

arima( y, order=c(0,0,2) )

## 
## Call:
## arima(x = y, order = c(0, 0, 2))
## 
## Coefficients:
##          ma1      ma2  intercept
##       0.7187  -0.1629    -0.0252
## s.e.  0.0310   0.0310     0.0541
## 
## sigma^2 estimated as 1.21:  log likelihood = -1514.87,  aic = 3037.74

ARMA(1,1)

zLogL=function(theta){
  n<-length(z)-1
  mu<-theta[1]
  phi<-theta[2]
  rho<-ARMAacf(ma=theta[3])
  gamma<-rho*theta[4]^2 * (1+theta[3]^2)
  A<-bandSparse(n,
               n,
               k=0:1,
               diagonals = cbind(rep(gamma[1],n),rep(gamma[2],n)),
               symmetric=TRUE)
  ii <- (1:n)+1
  logL<-log_mvnpdf(z[ii],mu +phi*z[ii-1],A)
}

optim( c(0,0,0,1), zLogL, method="BFGS" )

## $par
## [1] -0.01265133  0.94917268 -0.19805777  0.97966523
## 
## $value
## [1] 1397.103
## 
## $counts
## function gradient 
##       84       28 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

arima(z,order=c(1,0,1))

## 
## Call:
## arima(x = z, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.9497  -0.1974    -0.3693
## s.e.  0.0107   0.0335     0.4859
## 
## sigma^2 estimated as 0.9604:  log likelihood = -1399.69,  aic = 2807.38