# hostages.R
#
# Analysis of ITERATE monthly hostages data, 1968-2005
#
# Patrick T. Brandt


# Load the datafile
load("hostages-monthly.RData")

# Analyze the count properties and dynamics of the series -- source
# the PESTS code
source("../pests.r")


# Plot series and their ACFs -- this uses the counts.acf function from
# pests.r to compute the ACFs

# Also, add a vertical line to each plot for 9-11
par(mfrow=c(3,2))

plot(km[,"incident"], xlab="Date",
     ylab="", main="Monthly KIDNAP Events")
abline(v=as.yearmon(as.Date("2001-9-11")), lty=2, col="red")

plot(counts.acf(km[,"incident"], lag.max=12), main="ACF for Monthly KIDNAP Counts")


plot(sm[,"incident"], xlab="Date",
     ylab="", main="Monthly SKYJACK Events")
abline(v=as.yearmon(as.Date("2001-9-11")), lty=2, col="red")

plot(counts.acf(sm[,"incident"], lag.max=12), main="ACF for Monthly SKYJACK Counts")


plot(nkm[,"incident"], xlab="Date",
     ylab="", main="Monthly OTHER Events")
abline(v=as.yearmon(as.Date("2001-9-11")), lty=2, col="red")

plot(counts.acf(nkm[,"incident"], lag.max=12), main="ACF for Monthly OTHER Counts")


# Fit a simple PAR(p) model to each series

# Set the random number seed since the initial values for the
# estimators depend on random starting values.
set.seed(9860)

# Define lists of the datasets.
data.m <- c("km", "sm", "nkm")

# Set the maximum lag length for monthly analyses.

max.m <- 9

# Matrices to hold results
results.m <- matrix(0, length(data.m)*max.m, 4)
colnames(results.m) <- c("data", "p", "aic", "llf")

# Monthly analyses: loop over the datasets and the lag lengths
j <- 1

for (i in 1:length(data.m))
{
    for(p in 1:max.m)
    {
        # Set up smart starting values for the models
        if(p==1)
        { parp.init <- c(0.1) }
        else
        { if(inherits(tmp, "try-error")==FALSE)
              {parp.init <- c(tmp$coefs[1:(p-1),1], 0) }
          else
          { parp.init <- rep(0.1, p)}
      }
        # Estimate the PAR(p) model
        tmp <- try(Parp(get(data.m[i])[,1] ~ 1, p=p, parp.init=parp.init))

        # Assign the results into the results matrix and a new
        # object.  Note we have to test the return object here.

        assign(paste(data.m[i], p, sep="."), tmp)

        results.m[j,1] <- data.m[i]
        results.m[j,2] <- p

        if(inherits(tmp, "try-error")==FALSE)
        {
            results.m[j,3] <- tmp$aic
            results.m[j,4] <- tmp$llf
        }
        j <- j+1
    }
}



print(results.m)

# Now estimate the final models with the covariates included
# The models in the paper are more refined than these, so these are
# meant to be illustrative.

kidnaps.model <- Parp(km[,"incident"] ~ km[,"negotiationsuccess"] +
                      km[,"violentend"] + km[,"deathsincident"], p=8)

skyjacks.model <- Parp(sm[,"incident"] ~ sm[,"negotiationsuccess"] +
                       sm[,"violentend"] + sm[,"deathsincident"], p=5)

other.model <- Parp(nkm[,"incident"] ~ nkm[,"negotiationsuccess"] +
                    nkm[,"violentend"] + nkm[,"deathsincident"], p=4)

# Print out the results
print(kidnaps.model)
print(skyjacks.model)
print(other.model)

# Estimation of the dynamic multipliers for one unit changes in each
# covariate.  These compute the effects of changes in each variable
# holding the others at their means.

kidnaps.mults <- montecarlo.parp.multipliers(kidnaps.model,
                                             xvec=rep(1,4), n=1000)

skyjacks.mults <- montecarlo.parp.multipliers(skyjacks.model,
                                             xvec=rep(1,4), n=1000)

other.mults <- montecarlo.parp.multipliers(other.model, xvec=rep(1,4),
                                             n=1000)

# To summarize the simulated multipliers, you need to sum over the
# simulated values.  Here we summarize the 0.16, 0.5 and 0.84
# quantiles.  The columns of this output correspond to each of the
# covariates (with the intercept as the first).  The first matrix are
# the short run or impact multipliers, the second are the long run or
# total multipliers

apply(kidnaps.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))
apply(skyjacks.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))
apply(other.mults, c(2,3), quantile, probs=c(0.16, 0.5, 0.84))