# pests-example.r
#
# A simple example file for the PESTS source file.
#
# Patrick Brandt
#
# pbrandt@utdallas.edu
# School of Social Sciences
# University of Texas, Dallas
#
# This file contains no warranty or guarantee of performance.
# Use at your own risk for non commercial use
#

source("pests.r")

# Simulate some data.
#
#
# Generic parameters for the simulation....
#
set.seed(7680)
n <- 100                                  # Number of observations
x <- sqrt(0.5)*matrix(rnorm(n),ncol=1)    # Generic matrix of regressors

# PEWMA DGP parameters
#
# Prior mean for DGP is a0/b0 
a0 <- 3            # Prior for a_t
b0 <- 0.25          # Prior for b_t
w <- 0.7           # Omega, the EWMA parameter
d.pewma <- c(0.5)  # Delta, coefficent(s) on regressors

# generate data with no drift
y.pewma.nodrift <- pewmadgp(n,a0,b0,w,x,d.pewma)

# generate data with drift
x1 <- cbind(matrix(1, n, 1), x)
d.pewma.drift <- c(0.75, 0.5)

y.pewma.drift <- pewmadgp(n,a0,b0,w,x1,d.pewma.drift)

# PAR(1) DGP parameters
#
rho <- 0.5           # AR(1) coefficient
d.par1 <- c(1,0.5)   # Regressors (constant and coefficient)
y.par1 <- parpdgp(n,rho,x=matrix(cbind(rep(1,100),x)),d.par1)


# Plot the results of the dgps
# Note the trick to get Barlett corrected std errors on the ACFs

par(mfrow=c(3,2))

plot(y.pewma.nodrift, ylab="Count", main="PEWMA data without drift")
acftmp <- counts.acf(y.pewma.nodrift)
plot(acftmp)

plot(y.pewma.drift, ylab="Count", main="PEWMA data with drift")
acftmp <- counts.acf(y.pewma.drift)
plot(acftmp)

plot(y.par1, ylab="Count", main="PAR(1) data")
acftmp <- counts.acf(y.par1)
plot(acftmp)



# Run some models for comparisons.....
# Use the mnemonic of "dgp.estimator" to keep straight what we are
# doing....

pewma.Poisson <- glm(y.pewma.nodrift ~ 1 + x, family=poisson)
summary(pewma.Poisson)
par1.Poisson <- glm(y.par1 ~ x, family=poisson)
summary(par1.Poisson)

# Now use the Poisson results to start the PESTS models

# Notice that the drift term (intercept) is zero.
pewma.PEWMA <- Pewma(y.pewma.nodrift ~ x )

par1.PAR1 <- Parp(y.par1 ~ x)


# Now run some models with and without intercepts just as demos

Pewma.nodrift <- Pewma(y.pewma.nodrift ~ -1)
Pewma.nodrift2 <- Pewma(y.pewma.drift ~ -1, omega.init=0.2)
Pewma.drift <- Pewma(y.pewma.drift ~ 1, omega.init=Pewma.nodrift$coefs[1,1])