# # PURPOSE: Function for Poisson AR(1) analysis
# #
# # AUTHORS: Stefanos Kechagias, James Livsey, Vladas Pipiras
# #
# # DATE: April 2020
# #
# # VERSION: 3.6.3
#
#
# #---------simulate Poisson series---------#
# sim_poisson = function(n, ARMAmodel, lam){
# #====================================================================================#
# # PURPOSE Simulate Poisson with ARMA structure. See relation (1)
# # in https://arxiv.org/pdf/1811.00203.pdf
# #
# # INPUT
# # n series length
# # ARMAmodel list with ARMA parameters
# # lam marginal param
# # Output
# # X Poisson series
# #
# # Authors Stefanos Kechagias, James Livsey, Vladas Pipiras
# # Date July 2020
# # Version 3.6.3
# #====================================================================================#
#
# z = arima.sim(model = list(ar=ARMAmodel[[1]], ma=ARMAmodel[[2]]), n = n); z = z/sd(z)
# # The third argument of qnbinom is the prob of failure
# x = qpois(pnorm(z), lam)
# return(x)
# }
#
#
# sim_poisson_2 = function(n, ARMAmodel, lam){
# #====================================================================================#
# # PURPOSE Simulate Poisson with ARMA structure. See relation (1)
# # in https://arxiv.org/pdf/1811.00203.pdf Same as sim_poisson but also returns
# # the z
# #
# # INPUT
# # n series length
# # ARMAmodel list with ARMA parameters
# # lam marginal param
# # Output
# # X Poisson series
# #
# # Authors Stefanos Kechagias, James Livsey, Vladas Pipiras
# # Date July 2020
# # Version 3.6.3
# #====================================================================================#
# l = list()
# z = arima.sim(model = list(ar=ARMAmodel[[1]], ma=ARMAmodel[[2]]), n = n); z = z/sd(z)
# # The third argument of qnbinom is the prob of failure
# x = qpois(pnorm(z), lam)
#
# l[[1]] = z
# l[[2]] = x
# return(l)
# }
#
#
#
# # Generate AR series
# sim_pois_ar = function(n, phi, lam){
# #------------------------------------------------------------------------#
# # PURPOSE Simulate Poisson series with AR structure. See relation (1)
# # in https://arxiv.org/pdf/1811.00203.pdf
# #
# # INPUT
# # n series length
# # phi AR parameter
# # lam Marginal Parameter
# #
# # Output
# # x Poisson series
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #-----------------------------------------------------------------------#
#
# z = arima.sim(model = list(ar=phi), n = n); z = z/sd(z) # standardized
# x = qpois(pnorm(z), lam)
# return(x)
# }
#
#
#
#
#
#
# # # compute arbitrary polynomial
# # evalPolynomial_scalarX <- function(coefs, x){
# # #######################################################################
# # # PURPOSE compute arbitrary polynomial given coefficients and input
# # #
# # # INPUT
# # # coefs vector of coefficients
# # # x input
# # #
# # # Output
# # # value returns f(x) = coefs[1] * x^0 +
# # # coefs[2] * x^1 +
# # # ... +
# # # coefs[n] * x^(n-1)
# # #
# # # Authors Stefanos Kechagias, James Livsey
# # # Date January 2020
# # # Version 3.6.1
# # #######################################################################
# # if(length(coefs) == 1){
# # return(coefs) # handle scalar case
# # }
# # n <- length(coefs) - 1
# # out <- sum(coefs * x^(0:n))
# # return(out)
# # }
# #
# # # Evaluate k^th Hermite Polynomial at scalar x
# # evalHermPolynomial <- function(k, x){
# # #######################################################################
# # # PURPOSE Evaluate k^th Hermite Polynomial at scaler input.
# # # See relation (7) https://arxiv.org/pdf/1811.00203.pdf
# # #
# # # INPUT
# # # k index of Hermite coefficient
# # # x input to Hermite Polynomial
# # #
# # #
# # # Output
# # # value returns H_k(x). See relation (7)
# # # in https://arxiv.org/pdf/1811.00203.pdf
# # #
# # # Authors Stefanos Kechagias, James Livsey
# # # Date January 2020
# # # Version 3.6.1
# # #######################################################################
# # if(k < 0){
# # return(0)
# # }
# # coefs <- HermPolyCoefs[[k+1]]
# # out <- evalPolynomial(coefs, x)
# # return(out)
# # }
# #
# # # Evaluate Hermite Polynomial at a vector x
# # evalPolynomial <- Vectorize(evalPolynomial_scalarX, vectorize.args = "x")
#
# # # kth hermitte coefficient
# HermCoefP_k <- function(lam, k){
# #========================================================================#
# # PURPOSE Compute kth Hermite Coefficient. See relation (21) in
# # https://arxiv.org/pdf/1811.00203.pdf
# #
# # INPUT
# # lam Marginal parameter
# # k index of Hermite coefficient
# #
# # Output
# # HC_k kth hermite coeficient
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #========================================================================#
#
# # function for (k-1)st Hermite Polynomial
# her <- function(x){
# evalHermPolynomial(k-1, x)
# }
#
# # truncation numbe: check me
# N <- which(round(ppois(1:1000, lam), 7) == 1)[1]
#
# # compute terms in the sum of relation (21) in
# terms <- exp((-qnorm(ppois(0:N, lam, lower.tail= TRUE))^2)/2) *
# her(qnorm(ppois(0:N, lam, lower.tail = TRUE)))
#
# # take the sum of all terms
# HC_k <- sum(terms) / (sqrt(2*pi) * factorial(k))
# return(HC_k)
# }
#
# # # all hermitte coefficients
# HermCoefP <- function(lam){
# #========================================================================#
# # PURPOSE Compute all Hermite Coefficients. See relation (21) in
# # https://arxiv.org/pdf/1811.00203.pdf
# #
# # INPUT
# # lam Marginal parameter
# # maxCoef number of coefficients to return. Default = 20
# #
# # Output
# # HC All Hermite coeficients
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #========================================================================#
# h = 1:20 #check me
# HC = rep(NA, length(h)) # storage
# for(i in h) {
# HC[i] <- HermCoefP_k(lam = lam, k = i)
# }
#
# return(HC)
#
# }
#
# # Poisson AR(1) covariance matrix
# CovarPoissonAR = function(n,lam,phi){
# #========================================================================#
# # PURPOSE Compute the covariance matrix of a Poisson AR series.
# #
# # INPUT
# # lam Marginal parameter
# # phi AR parameter
# # n size of the matrix
# #
# # Output
# # GAMMA covariance matrix ofcount series
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #=========================================================================#
#
# # Hermite coeficients--relation (21) in https://arxiv.org/pdf/1811.00203.pdf
# HC = HermCoefP(lam)
#
# # ARMA autocorrelation function
# ar.acf <- ARMAacf(ar = phi, lag.max = n)
#
# # Autocovariance of count series--relation (9) in https://arxiv.org/pdf/1811.00203.pdf
# gamma_x = CountACVF(h = 0:(n-1), myacf = ar.acf, g = HC)
#
# # Final toeplitz covariance matrix--relation (56) in https://arxiv.org/pdf/1811.00203.pdf
# GAMMA = toeplitz(gamma_x)
# return(GAMMA)
# }
#
# # Poisson AR(1) Gaussian Lik
# GaussLogLik = function(theta, data){
# #-------------------------------------------------------------------------#
# # PURPOSE Compute Gaussian log-likelihood for Poisson AR series
# #
# # INPUT
# # theta parameter vector containing the marginal and AR parameters
# # data count series
# #
# # Output
# # loglik Gaussian log-likelihood
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #-------------------------------------------------------------------------#
#
# # retrieve parameters and sample size
# lam = theta[1]
# phi = theta[-1]
# n = length(data)
#
# # assign large likelihood value if not causal and if lambda outside range
# if(any(abs( polyroot(c(1, -phi)) ) < 1) || (lam < 0 || lam > 100) ){
# return(NA) #check me
# }
#
# # Compute the covariance matrix--relation (56) in https://arxiv.org/pdf/1811.00203.pdf
# GAMMA = CovarPoissonAR(n, lam, phi)
#
# # Compute the logdet and the quadratic part
# logLikComponents = EvalInvQuadForm(GAMMA, as.numeric(data), lam)
#
# # final loglikelihood value
# out = 0.5*logLikComponents[1] + 0.5*logLikComponents[2]
#
# # the following will match the above if you subtract N/2*log(2*pi) and don't multiply with 2
# # out = -2*dmvnorm(as.numeric(data), rep(lam, n), GAMMA, log = TRUE)
# return(out)
# }
#
# # wrapper to fit Poisson AR(1) using Gaussian Lik
# FitGaussianLik = function(initialParam, x){
# #######################################################################
# # PURPOSE Fit the Gaussian log-likelihood for Poisson AR series
# #
# # INPUT
# # initialParam parameter vector containing the marginal and AR parameters
# # x count series
# #
# # Output
# # optim.output$par parameter estimates
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #######################################################################
# optim.output <- optim(par = initialParam,
# fn = GaussLogLik,
# data = x,
# method = "BFGS",
# hessian=TRUE)
#
# nparms = length(initialParam)
# ParmEst = matrix(0,nrow=1,ncol=nparms)
# se = matrix(NA,nrow=1,ncol=nparms)
# loglik = rep(0,1)
#
# # save estimates, loglik and standard errors
# ParmEst[,1:nparms] = optim.output$par
# loglik = optim.output$value
# se[,1:nparms] = sqrt(abs(diag(solve(optim.output$hessian))))
#
# All = cbind(ParmEst, se, loglik)
# return(All)
#
# }
#
#
#
#
#
# GaussLogLik_fixedPhi = function(lam, phi_fixed, data){
# #-------------------------------------------------------------------------#
# # PURPOSE Compute Gaussian log-likelihood for Poisson AR series with Phi fixed
# #
# # INPUT
# # lam Poisson mean parameter
# # phi_fixed value to fix phi at
# # data count series
# #
# # Output
# # loglik Gaussian log-likelihood
# #
# # Authors Stefanos Kechagias, James Livsey
# # Date January 2020
# # Version 3.6.1
# #-------------------------------------------------------------------------#
#
# # retrieve parameters and sample size
# phi = phi_fixed
# n = length(data)
#
# # assign large likelihood value if not causal and if lambda outside range
# if(any(abs( polyroot(c(1, -phi)) ) < 1) || (lam < 0 || lam > 100)){
# return(NA) #check me
# }
#
# # Compute the covariance matrix--relation (56) in https://arxiv.org/pdf/1811.00203.pdf
# GAMMA = CovarPoissonAR(n, lam, phi)
#
# # Compute the logdet and the quadratic part
# logLikComponents = EvalInvQuadForm(GAMMA, as.numeric(data), lam)
#
# # final loglikelihood value
# out = 0.5*logLikComponents[1] + 0.5*logLikComponents[2]
#
# # the following will match the above if you subtract N/2*log(2*pi) and don't multiply with 2
# # out = -2*dmvnorm(as.numeric(data), rep(lam, n), GAMMA, log = TRUE)
# return(out)
# }
#
#
#
#
#
#
# # # count acvf for each h
# # CountACVF_h = function(h, myacf, g){
# # #=======================================================================#
# # # PURPOSE Compute the autocovariance matrix of the count series.
# # # See relation (9) in https://arxiv.org/pdf/1811.00203.pdf
# # # INPUT
# # # h acvf lag
# # # myacf autocorrelation of Gaussian series: rho_z
# # # g Hermitte Coefficients
# # # Output
# # # gamma_x count acvf
# # #
# # # Notes: See also CountACVF function--vectorized version of CountACVF_h
# # #
# # # Authors Stefanos Kechagias, James Livsey
# # # Date January 2020
# # # Version 3.6.1
# # #=======================================================================#
# #
# # k = length(g) #check me
# # gamma_x = sum(g^2 * factorial(1:k) * (myacf[h+1])^(1:k))
# # return(gamma_x)
# #
# # }
# #
# # # vectorized count acvf
# # CountACVF <- Vectorize(CountACVF_h, vectorize.args = "h")
# #
# # # fnc to evaluate Gaussial lik components
# # EvalInvQuadForm = function(A, v, DataMean ){
# # #=======================================================================#
# # # Evaluate quadrative form v`*inv(A)*v where
# # # A is symmetric positive definite
# # # Want QuadForm = v` * inv(A) * v = v` * w where A*w = v
# # # ==> (U`U)*w = v Cholesky decomp of A
# # # ==> First solve U` z = v for z,
# # # then solve U w = z for w */
# # #=======================================================================#
# #
# # U = chol(A)
# # z = forwardsolve(t(U), v - DataMean)
# # w = backsolve(U,z)
# # QuadForm = t(v-DataMean)%*%w
# #
# # logdetA = 2*sum(log(diag(U)))
# #
# # logLikComponents = c(logdetA, QuadForm)
# # return(logLikComponents)
# # }
#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.