#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

if( !require('rjags') ){
    install.packages('rjags', dependencies=TRUE,  repos='http://cran.rstudio.com/')
}
library(rjags)


save_plots = F

set.seed(0) 

library(R2WinBUGS)
library(BRugs)

# P 1: EPage 
#
data(CRSPmon,package="Ecdat")
ibm = CRSPmon[,2]
y = as.numeric(ibm)
N = length(y)
ibm_data = list("y","N")

inits = function(){
    list(mu=rnorm(1,0,0.3),tau=runif(1,1,10),nu=runif(1,1,30))
}

univt.mcmc = bugs(ibm_data,inits,
    model.file="Tbrate_t.bug",
    parameters=c("mu","tau","nu","sigma"),
    n.chains = 3, n.iter=2600, n.burnin=100,
    n.thin=1,
    program="OpenBUGS", 
    codaPkg=F, bugs.seed=1)

print(univt.mcmc,digits=4)
plot(univt.mcmc)


# Problem 2:
# 
mu = matrix(univt.mcmc$sims.array[,,1],ncol=1)
tau = matrix(univt.mcmc$sims.array[,,2],ncol=1)
nu = matrix(univt.mcmc$sims.array[,,3],ncol=1)
sigma = matrix(univt.mcmc$sims.array[,,4],ncol=1)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_2_ts_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,2))
ts.plot(mu,xlab="iteration", ylab="", main="mu") 
ts.plot(sigma,xlab="iteration", ylab="", main="sigma")
ts.plot(nu,xlab="iteration", ylab="", main="df") 
if( save_plots ){ dev.off() }

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_2_acf_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,2))
acf(mu,main="mu")
acf(sigma,main="sigma")
acf(nu,main="df")
if( save_plots ){ dev.off() }

library("fGarch") # for skewness and kurtosis
skewness(nu,method="moment") 
kurtosis(nu,method="moment")

# Lets check that for a t distribution with known df the kurtosis using the above function equals 3*(nu-2)/(nu-4):
df = 16
samps = rt( 500000, df=df )
c( kurtosis(samps, method="moment"), 3*(df-2)/(df-4) ) 

# Problem 3
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_2_hist_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,2))
hist(mu,main="mu")
hist(sigma,main="sigma")
hist(nu,main="df")
if( save_plots ){ dev.off() }


if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_2_density_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,2))
plot(density(mu),main="mu")
plot(density(sigma),main="sigma")
plot(density(nu),main="df")
if( save_plots ){ dev.off() }

c( skewness(mu), skewness(sigma), skewness(nu) )


# Problem 4
# 
nu_too_small = nu <= 4 # check for nu's that give infinite kurtosis 
posterior_kurtosis = 3*(nu-2)/(nu-4)
posterior_kurtosis[nu_too_small] = +Inf

quantile( posterior_kurtosis, c( 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 ) ) 


library(boot)
boot.fn = function(data,index){
    return( kurtosis(data[index],method="moment") )
}
boot.kurtosis = boot( y, boot.fn, R=1000 ) 

quantile( boot.kurtosis$t, c( 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 ) ) 

boot.ci( boot.kurtosis, type="perc", conf=0.9 ) 


# Problem 5:
# 
TbGdpPi = read.csv('../../BookCode/Data/TbGdpPi.csv', header=TRUE)
# r = the 91-day Treasury bill rate
# y = the log of real GDP
# pi = the inflation rate
#
TbGdpPi = ts(TbGdpPi, start=1955, freq=4)
Tbrate = TbGdpPi
del_dat = diff(Tbrate)
y = as.numeric(del_dat[,2])
y = y - mean(y)
N = length(y)
GDP_data = list("y","N")

#### AR 1, GDP data #####
inits = function(){ list(phi=rnorm(1,0,0.3), tau=runif(1,1,10)) }
ar1.mcmc = bugs(GDP_data, inits,
    model.file="ar1.bug",
    parameters=c("phi", "sigma"),
    n.chains=3, n.iter=2600, n.burnin=100,
    n.thin=1,
    program="OpenBUGS",
    bugs.seed=1)

print(ar1.mcmc,digits=3)
plot(ar1.mcmc)
arima(y,order=c(1,0,0))

# Plot time series of the parameters phi and sigma:
#
phi = matrix(ar1.mcmc$sims.array[,,1], ncol=1)
sigma = matrix(ar1.mcmc$sims.array[,,2], ncol=1)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_5_phi_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
ts.plot(phi, main="phi")
acf(phi, main="phi ACF")
if( save_plots ){ dev.off() }

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter20/chap_20_prob_5_sigma_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
ts.plot(sigma, main="sigma")
acf(sigma, main="sigma ACF")
if( save_plots ){ dev.off() }



if( FALSE ){
    # "not finished" 
    stopifnot( FALSE )
    
    #### MA 1, simulated data #####
    library(R2WinBUGS)

    set.seed(5640)
    N = 600
    y = arima.sim(n=N, list(ma=-0.5), sd=0.4)
    y = as.numeric(y)
    q = 5
    ma.sim_data = list("y","N","q")


    inits.ma = function(){ list(theta=rnorm(1,-0.5,0.1), tau=runif(1,5,8)) }
    ma1.mcmc = bugs( ma.sim_data, inits.ma,
        model.file="ma1.bug",
        parameters=c("theta", "sigma", "ypred"),
        n.chains=3, n.iter=3000, n.burnin=500,
        n.thin=1,
        program="OpenBUGS",
        codaPkg=F, bugs.seed=1)

    print(ma1.mcmc, digits=3)
    plot(ma1.mcmc)

    # Problem 7
    #
    set.seed(5640)
    N = 600
    y = arima.sim( n=N, list(ar = 0.9, ma=-0.5), sd=0.4 )
    y = as.numeric(y)

}