#-- If not installed on PC ---
#  install.packages("tseries")
#  install.packages("xts")
#  install.packages("forecast")
#


#-- To load this page --
#
#    source("https://nmimoto.github.io/R/TS-00.txt")
#
#  Above source() automatically installs:
#    Randomness.tests()
#    Stationarity.tests()
#    Rolling1step.forecast()


library(xts)
library(forecast)


#---------------------------------------------
Randomness.tests <- function( A, plott=TRUE) {


  # remove NA
  if (sum(is.na(A))>0) A <- A[which(is.na(A)==0)] 

  # turn to numeric if xts
  if (is.xts(A)) A <- as.numeric(A)  
  
  library(tseries) 

  L1 <- Box.test(A,   lag = 15, type = "Ljung-Box")
  L2 <- Box.test(A,   lag = 20, type = "Ljung-Box")
  L3 <- Box.test(A,   lag = 25, type = "Ljung-Box")
  L4 <- Box.test(A^2, lag = 15, type = "Ljung-Box")
  L5 <- Box.test(A^2, lag = 20, type = "Ljung-Box")
  L7 <- jarque.bera.test(A)
  S1 <- sd(A)

  if (plott) {

    layout( matrix(c(rep(1,4),2,8,3,4,5,6,7,8), 2, 6, byrow=T) )

    par(oma = c(2,2,4,2) + 0.1, mar = c(4,1,1,1) + 0.1 )

    plot(A, type='l', xlab="Series")

    plot(density(A, bw="SJ-ste"), main= "", xlab="KDE")

    acf( A, main="", xlab="ACF")
    pacf(A, main="", xlab="PACF")
    acf(abs(A), main="", xlab="ACF of abs(A)")
    acf(A^2,    main="", xlab="ACF of A^2")
    qqnorm(A,    main="", xlab="Normal q-q")
 
    
    plot(c(-1,-1), xlim=c(0,1), ylim=c(0,1), ann=F, axes=F)

    text( 0.5,0.96, paste("Box-Ljung test"), cex=1.3 )
    text( 0.5,0.9,  paste("H=15:p=", round(L1$p.value, 3)), cex=1.2 )
    text( 0.5,0.85, paste("H=20:p=", round(L2$p.value, 3)), cex=1.2 )
    text( 0.5,0.8,  paste("H=25:p=", round(L3$p.value, 3)), cex=1.2 )
    text( 0.5,0.71, paste("McLeod-Li test") , cex=1.3 )
    text( 0.5,0.65, paste("H=15:p=", round(L4$p.value, 3)), cex=1.2 )
    text( 0.5,0.6,  paste("H=20:p=", round(L5$p.value, 3)), cex=1.2 )
    text( 0.5,0.35, paste("Jaque-Bera test") , cex=1.3 )
    names(L7$p.value) <- ""
    text( 0.5,0.3,  paste("p= ", round(L7$p.value, 4)), cex=1.2 )
    text( 0.5,0.2,  paste("SD ", round(S1, 4)), cex=1.2 )

    mtext(side = 3, line=1, outer=T, text = "Randomness.tests()", cex=1.5) #-- add big title to panel plot
    
    layout( matrix(1, 1, 1) )
  }

  P <- matrix(c(round(L1$p.value, 3),round(L2$p.value,3),round(L3$p.value,3),round(L4$p.value,3),
                round(L5$p.value,3), round(L7$p.value,3),round(S1, 3) ), 1, 7)
 
  colnames(P) <- c("BL15","BL20","BL25","ML15","ML20","JB","SD") 

  cat("  B-L test H0: the series is uncorrelated\n")
  cat("  M-L test H0: the square of the series is uncorrelated\n")
  cat("  J-B test H0: the series came from Normal distribution\n")
  cat("  SD         : Standard Deviation of the series\n\n")

  return(P)
}




#---------------------------------------------
Stationarity.tests <- function( A ) {

  library(tseries) 

  if (sum(is.na(A))>0){  A <- A[which(is.na(A)==0)]  }

  L1 <- adf.test(A)
  L2 <- pp.test(A)
  L3 <- kpss.test(A)

  P <- matrix(c(round(L3$p.value,3),round(L1$p.value, 3),round(L2$p.value,3)), 1, 3)
  colnames(P) <- c("KPSS","ADF","PP")
  rownames(P) <- c("p-val:")

  return(P)
}



#--- Usage Example
#
#  x = rnorm(200)
#  Randomness.tests(x)
#  Stationarisy.tests(x)
#




#---------------------------------------------
MV.av <- function(X,m){
  # Moving avegrate that uses past m observations.  Only takes xts object.
  # Output $mean and $sd, and they are also xts.
  # If m=20, 1st $mean is mean(X[1:20]), which is indexed at 21st day.
  
    if (is.xts(X)) {
    
        n <- length(X)
        
        Mav1 <- xts(matrix(0,n-m,1), order.by=index(X[(m+1):n]))
        Mav2 <- xts(matrix(0,n-m,1), order.by=index(X[(m+1):n]))
        
        for (i in  1:(n-m)){
            Mav1[i] <-  mean(X[(i-1)+1:m])
            Mav2[i] <-  sd(X[(i-1)+1:m])   
        }
        
    } else {
    
        print("X needs to be xts")
        
    }
    
    return(list("mean"=Mav1, "sd"=Mav2))
    
}
    
    



#--- Normality Test ---------------------------
KStest <- function(X) {

    n <- length(X)
    X <- sort(X)
    F1 <- (1:n)/n
    F2 <- (0:(n-1))/n
    a <- max(abs( F1 - F0(X) ))
    b <- max(abs( F2 - F0(X) ))
    c <- max(a,b)*sqrt(n)
    return(c)
    
}


#-- built-in ks.test() lacks the sqrt(n) scaling!!
#
#    ks.test(X, "pnorm", 0, 1)$statistic*sqrt(n)
#




#--- Rolling 1-step version 2022.3 ---------------------------
Rolling1step.forecast <- function(Y,
                                  window.size,
                                  Arima.order,
                                  include.mean=TRUE,
                                  include.drift=FALSE,
                                  lambda=NULL,           # NULL=no transformaton. 0=Log
                                  xreg=NULL,             # NULL=regression.  TRUE=Linear Trend is present
                                  seasonal=c(0,0,0),     # seasonal component
                                  show_plot=TRUE){
  
  
  library(forecast)
  # take xreg=FALSE/NULL as xreg=NULL
  if (is.null(xreg)==1 || xreg==0){ xreg = NULL }
  
  Yhat  <- Yhat.CIu <- Yhat.CIl <- 0        # initialize what needs to be saved
  for (i in 1:(length(Y)-window.size)) {
    
    
    if (is.null(xreg)){ xreg=NULL } else { xreg=time(Y[i:(i+window.size-1)]) }
    
    # Use Arima() with given options
    Fit00 <- tryCatch(
      Arima( Y[i:(i+window.size-1)],
             order=Arima.order,
             include.mean = include.mean,
             include.drift = include.drift,
             lambda = lambda,
             seasonal = seasonal, 
             xreg=xreg), 
      error = function(ss){
        cat(paste("    i=",i,"  MLE-CSS failed.  Using CSS.\n"))
        return(
          Arima( Y[i:(i+window.size-1)],
                 order=Arima.order,
                 include.mean = include.mean,
                 include.drift = include.drift,
                 lambda = lambda,
                 seasonal = seasonal,
                 xreg=xreg,
                 method="CSS") ) } )
    
    # Model one step forecast  
    if (is.null(xreg)){ 
      Y.h  <- forecast(Fit00, 1)
    } else {
      Y.h  <- forecast(Fit00, 1, xreg=last(time(Y[i:(i+window.size-1)]))+1/frequency(Y))
    }
    Yhat[i]     <- Y.h$mean
    Yhat.CIu[i] <- Y.h$lower[2]
    Yhat.CIl[i] <- Y.h$upper[2]     
  }
  
  # Yhat starts at window.size+1 up to length(Y)
  Yhat     = ts(Yhat,     start=time(Y)[window.size+1], freq=frequency(Y))
  Yhat.CIu = ts(Yhat.CIu, start=time(Y)[window.size+1], freq=frequency(Y))
  Yhat.CIl = ts(Yhat.CIl, start=time(Y)[window.size+1], freq=frequency(Y))
  Y1 = window(Y, end=time(Y)[window.size],     freq=frequency(Y))
  Y2 = window(Y, start=time(Y)[window.size+1], freq=frequency(Y))
  
  #- Calculate prediction performance
  Pred.error = Y2 - Yhat
  Pred.rMSE =  sqrt(  mean( (Pred.error)^2 ) )     # prediction root Mean Squared Error
  
  cat(paste("\n  Total length",length(Y),", window size",window.size,".\n  Last", 
            (length(Y)-window.size), "obs retrospectively forecasted with Rolling 1-step\n", 
            "     prediction using same order and fized window size.\n"))
  cat(paste("\n  Average Prediction Error: ", round(mean(Pred.error), 4)))
  cat(paste("\n  root Mean Squared Error of Prediction:  ", round(Pred.rMSE, 4)), "\n\n")
  
  if (show_plot) {
    #- Plot the prediction result with original data
    layout(matrix(c(1,1,1,1,1,1,2,3,4), 3, 3, byrow=TRUE))    # 3 plots at once
    
    plot(  Y,        type="n", ylab="", main="Rolling 1-step prediction (Red=Actual, Blue=Prediction)")  
    lines( Y1,       type="o") 
    lines( Y2,       type="o", col="red")
    lines( Yhat,     type="o", col="blue")
    lines( Yhat.CIu, type="l", col="gray30", lty=2)
    lines( Yhat.CIl, type="l", col="gray30", lty=2)
    
    # - Plot the prediction error
    plot( (window.size+1):length(Y), Pred.error, type="o", xlab="", main="Prediction Error (Red-Blue)" ) 
    abline(h=0, lty=2)
    
    acf(Pred.error, main="ACF of Pred.error")
    layout(1)                    # reset the layout
  }
  
  P <- matrix(c(round(mean(Pred.error), 4), round(Pred.rMSE, 4)), 1, 2)
  colnames(P) <- c("mean pred error","rMSE") 
  output = list(Actual=Y2, Predicted=Yhat)
  return(output)
  
  # Rolling 1-step version 2022.3
}


# Example of how to use it
# Pred = Rolling1step.forecast(Y, window.size=50, Arima.order=c(1,0,1), include.mean=TRUE)
# Pred
#
####################