Nothing
## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
cond_loglik <- function(model, x, index){ # todo: the same for the derivative?
if(missing(index)) # or for the derivative of the
index <- (max(model@order) + 1) : length(x) # one in the M-step of EM.
wrk <- mix_pdf(model, x, index)
wrk <- log(wrk)
wrk <- sum(wrk)
wrk
}
## For use with seasonal coefficients.
## todo: inefficient, could adapt mix_pdf etc.. to make it quicker
cond_loglikS <- function(model, x, index){
n <- length(x)
g <- length(model@prob)
if(missing(index))
index <- (max(model@arcoef@ps * model@arcoef@s)+1):length(x)
yhat <- matrix(0, nrow = length(index) , ncol = g) # rows contain same time,
# columns contain same component
for(t in seq_along(index)){ ## Calculate predictions
for(k in 1:g){
yhat[t,k] <- model@shift[k] +
sum(model@arcoef@a[[k]] * x[(index[t] - 1):(index[t] - model@order[k])]) +
sum( model@arcoef@as[[k]] * ui(x,index[t], model@arcoef@s *
(1:model@arcoef@ps[k])) )
}
}
pdf <- numeric(length(index)) ## Calculate density
for(t in seq_along(index)){
for(k in 1:g){
pdf[t] <- pdf[t] + model@prob[k] * dnorm(x[index[t]], yhat[t, k], model@scale[k])
}
}
ll <- log(pdf) ## loglikelihood
sum(ll)
}
lik_params_bounds <- function(model){
k <- .nmix(model) # number of mixture components
n <- length(lik_params(model)) # effective # of parameters
lower <- rep(-Inf, n) # lower bounds
lower[1:(k-1)] <- rep(0,k-1) # probabilities are >= 0
lower[k:(k+k-1)] <- rep(0,k) # st.dev's are >=0
upper <- rep(Inf, n)
upper[1:(k-1)] <- rep(1,k-1) # probabilities are <= 1
list(lower = lower, upper = upper)
}
## rdist - list of functions,
## z - vector of positive integers ('regimes')
## return value - numeric vector of same length as z
mixARnoise_sim <- function(rdist, z){
g <- length(rdist) # g must be positive
n <- length(z)
if(g == 1){
do.call(rdist[[1]], list(n)) # the dist is one and the same, z not needed in this case
}else{
wrk <- rep(as.numeric(NA), n)
for(i in 1:g){
ni <- sum(z == i)
wrk[z == i] <- do.call(rdist[[i]], list(ni))
}
wrk
}
}
## todo: tova e nay-dobre da stane generic.
## todo: change the default for nskip to 0?
## set default for init using the unconditional mean?
mixAR_sim <- function(model, n, init, nskip = 100, flag = FALSE){
sigma <- model@scale
shift <- model@shift
zdist <- model@prob
K <- length(zdist) # number of components
## 'seas_flag' is used below at places were the seasonal case needs its own treatment
seas_flag <- inherits(model@arcoef, "raggedCoefS")
if(seas_flag){
s <- model@arcoef@s
order <- matrix(0, nrow = K, ncol = max(model@order + model@arcoef@ps))
for(i in seq_len(K)){
order[i, 1:(model@order[i] + model@arcoef@ps[i])] <-
c(seq_len(model@order[i]), s * seq_len(model@arcoef@ps[i]))
}
}else
order <- model@order
if(length(init) < max(order))
stop("The length of argument 'init' must be at least equal to the maximal order.")
resinit <- tail(init, max(order)) # lastn(init, max(order))
# use only the last max(order) values in init.
ntot <- nskip + n
nresinit <- length(resinit)
z <- sample.int(K, size = ntot, replace = TRUE, prob = zdist)
noisedist <- noise_rand(model)
res <- mixARnoise_sim(noisedist, z)
res <- c(resinit, res) # now prepend the initial values
z <- c(numeric(nresinit), z)
ar <- model@arcoef
## needs efficient implementation, this is extremely inefficient (but straightforward)
for(i in nresinit + (1:ntot)){
k <- z[i]
wrk <- sigma[k] * res[i] + shift[k]
# 2018-12-09: added simulation from seasonal AR
if(seas_flag){
filt <- c(ar@a[[k]], ar@as[[k]])
for(j in seq_along(order[k, which(order[k, ] !=0)])){
wrk <- wrk + filt[j] * res[i - order[k, j]]
}
} else{
filt <- ar[[k]]
for(j in seq_len(order[k])){ # 2011-07-12: was 1:order[k]
wrk <- wrk + filt[j] * res[i - j]
}
}
res[i] <- wrk
}
if(flag) structure(lastn(res, n), regimes = lastn(z, n)) # return the last n values (could
else lastn(res, n) # use tail()); if 'flag', attach
} # 'z' as attribute to the result.
# 2012-12-03: new
# quick fix for Shahadat, consolidiray later.
mixAny_sim <- function(model, n, init, nskip=100, flag = FALSE,
theta, # ma coef, a list
galpha0, # alpha0[k], k=1,...,g.
galpha, # garch alpha
gbeta # garch beta
){
order <- model@order
sigma <- model@scale
shift <- model@shift
zdist <- model@prob
g <- length(model@prob)
thetaorder <- sapply(theta, length)
alphaorder <- sapply(galpha, length)
betaorder <- sapply(gbeta, length)
baseind <- max(betaorder) + 1 # baseind - j is tauinv[t-j]
# column k is for the k-th component
tauinvmat <- matrix(rep(galpha0), byrow = TRUE, ncol = g, nrow = max(betaorder))
K <- length(zdist) # number of components
if(length(init) < max(order))
stop("The length of argument 'init' must be at least equal to the maximal order.")
resinit <- tail(init, max(order,alphaorder,thetaorder)) # lastn(init, max(order))
# use only the last max(order) values in init.
ntot <- nskip + n
nresinit <- length(resinit)
z <- sample.int(K, size=ntot, replace=TRUE, prob=zdist)
noisedist <- noise_rand(model)
res <- mixARnoise_sim(noisedist, z)
eps <- c(numeric(nresinit), res) # needed for the moving averages
res <- c(resinit, res) # now prepend the initial values
z <- c(numeric(nresinit), z)
ar <- model@arcoef
## needs efficient implementation, this is extremely inefficient (but straightforward)
for(i in nresinit+(1:ntot)){
k <- z[i]
tauinv <- galpha0 # compute new tauinv # eq. (5.2), Shahadat
for(icomp in 1:g){
for(j in seq_len(alphaorder[icomp])){
tauinv[icomp] <- tauinv[icomp] + galpha[[icomp]][j]*eps[i-j]^2
}
for(j in seq_len(betaorder[icomp])){
tauinv[icomp] <- tauinv[icomp] +
gbeta[[icomp]][j]*tauinvmat[baseind - j, icomp]
}
}
sigma[k] <- sqrt(tauinv[k]) # multiply epsilon by this below
tauinvmat <- rbind(tauinvmat[-1, ], tauinv) # need only the last few tauinv's
# so drop the oldest
filt <- ar[[k]]
wrk <- sigma[k]*res[i] + shift[k]
for(j in seq_len(order[k])){ # 2011-07-12: was 1:order[k]
wrk <- wrk + filt[j]*res[i-j]
}
maco <- theta[[k]] # ma part
for(j in seq_len(length(maco))){
wrk <- wrk + maco[j]*eps[i-j]
}
res[i] <- wrk
}
if(flag) structure(lastn(res,n), regimes = lastn(z,n)) # return the last n values (could
else lastn(res,n) # use tail()); if 'flag', attach
} # 'z' as attribute to the result.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.