# straits.R
#
# Analysis of militarized interstate disputes (MIDS) for China and
# Taiwan, 1950-2001
#
# 20100627 : Initial version
#
#

# Load libraries and code
library(zoo)
source("../pests.r")

# Load the data
load("StraitsMIDS.RData")

# Make d into a ts() object
d <- as.ts(d)

# Plot the data series
pdf(file="CHNTWNMIDS.pdf", width=6, height=6)
plot(d,
     main="China and Taiwan threats and actions, MIDS: 1950-2001")
dev.off()

# Compute the ACFs for the standardized data

pdf(file="CHNTWNMIDS-acf.pdf", width=6, height=6)
acf(apply(d, 2, scale))
dev.off()


# Fit the Poisson and negative binomial models
CT.psn <- glm(d[,1] ~ 1, family=poisson())
CA.psn <- glm(d[,2] ~ 1, family=poisson())

TT.psn <- glm(d[,3] ~ 1, family=poisson())
TA.psn <- glm(d[,4] ~ 1, family=poisson())

library(MASS)
CT.nb <- glm.nb(d[,1] ~ 1)
CA.nb <- glm.nb(d[,2] ~ 1)

TT.nb <- glm.nb(d[,3] ~ 1)
TA.nb <- glm.nb(d[,4] ~ 1)

summary(CT.psn)
summary(CA.psn)
summary(TT.psn)
summary(TA.psn)

summary(CT.nb)
summary(CA.nb)
summary(TT.nb)
summary(TA.nb)


# Now fit the PAR(p) models
CT.par1 <- Parp(d[,1] ~ 1, p=1)
CT.par2 <- Parp(d[,1] ~ 1, p=2)
CT.par3 <- Parp(d[,1] ~ 1, p=3)

CA.par1 <- Parp(d[,2] ~ 1, p=1)
CA.par2 <- Parp(d[,2] ~ 1, p=2)

# Use the PAR(2) results to initialize the PAR(3) model.
CA.par3 <- Parp(d[,2] ~ 1, p=3, parp.init=c(0.6,0.6,0))

# Taiwan threats data are not dynamic -- so these models will not
# converge.
#TT.par1 <- Parp(d[,3] ~ 1, p=1, parp.init=0, init.param=-2)
#TT.par2 <- Parp(d[,3] ~ 1, p=2)
#TT.par3 <- Parp(d[,3] ~ 1, p=3)

TA.par1 <- Parp(d[,4] ~ 1, p=1)
TA.par2 <- Parp(d[,4] ~ 1, p=2)

# Now compute the dynamic multipliers using parp.impulse() over 10
# years.  Use the best fitting models based on AIC from above.

n <- 10
CT.irf <- parp.impulse(CT.par2, x.impulse=matrix(rep(0,n), ncol=1), n)
CA.irf <- parp.impulse(CA.par2, x.impulse=matrix(rep(0,n), ncol=1), n)
TA.irf <- parp.impulse(TA.par1, x.impulse=matrix(rep(0,n), ncol=1), n)

# Now plot the impulses
irfs <- cbind(CT.irf, CA.irf, TA.irf)
colnames(irfs) <- colnames(d)[c(1,2,4)]

pdf(file="impulses.pdf", width=6, height=6)
plot(ts(irfs), plot.type="single", lwd=2, col=1:3, ylab="Impulse Response")
abline(h=0)
legend(2, -2, colnames(irfs), lwd=2, col=1:3)
dev.off()