# HW 3 Q 6 JLCY, 19 Dec 2013

############################################################ 

### Simulate from the raw spectrum ############

# clear variables
rm(list = ls())

par(ask = F)

library(tseries)

## Warning: package 'tseries' was built under R version 3.0.2

# source('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/rawspec.R')
# source('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/smoothspec.R')

source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\rawspec.r")
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\smoothspec.r")

Qnumber = 12

source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\loaddata.r")

## Package `sm', version 2.2-5: type help(sm) for summary information

y = May

deltaT = 1
t = 1:length(y)

############################################################ 

##### Get the ak and bk and it is equivalent to to the FFT - see below

# number of frequencies possible..
N = length(y)
freqs = 2 * pi * (1:round((N/2)))/N
plotfreqs = freqs/(2 * pi * (deltaT))
nfreq = length(freqs)
a = 1:nfreq
b = 1:nfreq

###################### 

# 1. Fit a best AR model to the data
arbest = ar((y - mean(y)), order.max = 25)
p = order(arbest$aic)[1]

## AR-1 model
armodel = arima((y - mean(y)), order = c(p, 0, 0), include.mean = TRUE, method = "ML")

############################################################ 

# 2. Simulate from the fitted model of the same length

## matrix to hold the raw spectrum..
specsimrawar = matrix(0, ncol = nsim, nrow = nfreq)
specsimrawsur = matrix(0, ncol = nsim, nrow = nfreq)

### matrix to hold the smoothed spectrum
specsimrawars = matrix(0, ncol = nsim, nrow = nfreq)
specsimrawsurs = matrix(0, ncol = nsim, nrow = nfreq)

# May obs statistics
ymean = mean(y)
ysd = sd(y)

# statistics of AR model
armean = c(1:nsim)
arvar = c(1:nsim)
arcor = c(1:nsim)
arskw = c(1:nsim)
armax = c(1:nsim)
armin = c(1:nsim)

# spectrum model
simpdfs = matrix(0, nrow = nsim, ncol = nevals)

# statistics of AR model
armeans = c(1:nsim)
arvars = c(1:nsim)
arcors = c(1:nsim)
arskws = c(1:nsim)
armaxs = c(1:nsim)
armins = c(1:nsim)

for (isim in 1:nsim) {

    ## simulate from AR1 model
    ysim = arima.sim(n = N, list(ar = c(armodel$coef[1:p])), sd = sqrt(armodel$sigma2)) + 
        ymean

    # put the AR1 array into a matrix
    simdismon = matrix(ysim, nrow = 1)

    simdismon = t(simdismon)

    maysim = simdismon

    simpdf[isim, ] = sm.density(maysim, eval.points = xeval, display = "none")$estimate

    ############################################################ 

    ## simulate from spectrum

    ysims = surrogate(y, fft = TRUE, amplitude = T)
    # ysims = surrogate(y,fft=T,amplitude=F)

    simdismons = matrix(ysims, nrow = 1)

    simdismons = t(simdismons)

    maysims = simdismons

    simpdfs[isim, ] = sm.density(maysims, eval.points = xeval, display = "none")$estimate

    # a = density(maysims, n = 100)

    # simpdfs[isim,] = a$x

    # simpdfs[isim,] = sm.density(maysims, display='none')$estimate

    source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\sim_basicstatQ6.r")

    ############################################################ 

    # 3. Compute the raw spectrum of the simulations

    specsimrawar[, isim] = rawspec(ysim, deltaT)
    specsimrawsur[, isim] = rawspec(ysims, deltaT)

    ### Compute the smooth spectrum of the simulations..

    specsimrawars[, isim] = smoothspec(ysim, deltaT)
    specsimrawsurs[, isim] = smoothspec(ysims, deltaT)

}

############################################################ 

## spectrum of the original data

#### raw spectrum of original data
pgram = rawspec(y, deltaT)

#### smooth spectrum of the data
pgrams = smoothspec(y, deltaT)

### Plot the raw spectrum from AR simulations

plot(plotfreqs, pgram, xlab = "Freq. cy/yr", ylab = "Spectrum", ylim = range(specsimrawar, 
    pgram), col = "red", lwd = 2)

for (i in 1:nsim) {
    lines(plotfreqs, specsimrawar[, i], col = "grey")
}

# lines(plotfreqs,pgram,col='red',lwd=3)

############################################################ 

par(ask = TRUE)

## plot the raw spectrum from surrogate data - i..e, simulations from
## spectrum

plot(plotfreqs, pgram, xlab = "Freq. cy/yr", ylab = "Spectrum", ylim = range(specsimrawsur, 
    pgram), col = "red", lwd = 3, type = "l")

for (i in 1:nsim) {
    lines(plotfreqs, specsimrawsur[, i], col = "grey")
}

lines(plotfreqs, pgram, col = "red", lwd = 3)

############################################################ 

############## Plot the smooth spectrum of the AR simulations

plot(plotfreqs, pgrams, xlab = "Freq. cy/yr", ylab = "Spectrum", ylim = range(specsimrawars, 
    pgrams), col = "red", lwd = 3, type = "l")

for (i in 1:nsim) {
    lines(plotfreqs, specsimrawars[, i], col = "grey")
}

lines(plotfreqs, pgrams, col = "red", lwd = 3)

############################################################ 

############## Plot the smooth spectrum of the surrogate data

plot(plotfreqs, pgrams, xlab = "Freq. cy/yr", ylab = "Spectrum", ylim = range(specsimrawsurs, 
    pgrams), col = "red", lwd = 3, type = "l")

for (i in 1:nsim) {
    lines(plotfreqs, specsimrawsurs[, i], col = "grey")
}

lines(plotfreqs, pgrams, col = "red", lwd = 3)

############################################################ 

# plot may flow from AR model
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\mayflow.r")

# boxplots from AR model
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\boxplots6Q6.r")

simpdf = simpdfs

# plot may flow from spectrum
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\mayflow.r")

armean = armeans
arvar = arvars
arcor = arcors
arskw = arskws
armax = armaxs
armin = armins

# boxplots from spectrum
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\boxplots6Q6.r")

############################################################