# R/rma.R In freqdom: Frequency Domain Based Analysis: Dynamic PCA

#### Documented in rma

#' Generates a zero mean vector moving average process.
#'
#' This simulates a vector moving average process
#' \deqn{
#'   X_t=\varepsilon_t+\sum_{k \in lags} \Psi_k \varepsilon_{t-k},\quad 1\leq t\leq n.
#' }
#' The innovation process \eqn{\varepsilon_t} is either multivariate normal or multivarite \eqn{t} with
#' a predefined covariance/scale matrix sigma and zero mean. The noise is generated with the
#' package \code{mvtnorm}. For Gaussian noise we use \code{\link[mvtnorm]{rmvnorm}}. For Student-t noise we use
#' \code{\link[mvtnorm]{rmvt}}. The parameters \code{sigma} and \code{df} are imported as arguments, otherwise we use default settings.
#'
#' @title Moving average process
#' @param n number of observations to generate.
#' @param d dimension of the time series.
#' @param Psi a \code{\link{timedom}} object with operators \code{Psi$operators}, where \code{Psi$operators[,,k]}
#' is the operator on thelag \code{lags[k]}. If no value is set then we generate a vector moving average process
#' of order \eqn{1}. Then, \code{Psi$lags = c(1)} and \code{Psi$operators[,,1]} is proportional to \eqn{\exp(-(i+j)\colon 1\leq i, j\leq d)} and such
#' that the spectral radius of \code{Psi[,,1]} is \eqn{1/2}.
#' @param noise \code{mnormal} for multivariate normal noise or \code{mt} for multivariate \eqn{t} noise. If not specified \code{mnormal} is chosen.
#' @param sigma covariance  or scale matrix of the innovations. If NULL then the identity matrix is used.
#' @param df degrees of freedom if \code{noise = "mt"}.
#' @return A matrix with \code{d} columns and \code{n} rows. Each row corresponds to one time point.
#' @export
#' @keywords simulations
rma = function(n, d = 2, Psi = NULL, noise = c("mnormal","mt"), sigma = NULL, df = 4)
{
if (!is.null(Psi))
d = dim(Psi)[1]

if (n < 1)
stop ("n must be a positive integer")
if (d < 1)
stop ("d must be a positive integer")
if (length(noise)>1)
noise = 'mnormal'
if (is.null(sigma))
sigma = diag(d)

if (det(sigma)<0 || !isSymmetric(sigma))
stop("sigma is not a covariance matrix.")

# if no operator then make some default
if (is.null(Psi)){
Psi = exp(-(1:d))%*%t(exp(-(1:d)))
Psi = Psi/norm(Psi,type="2")/2
}

if(is.array(Psi)) nlags = dim(Psi)[3]
if(is.matrix(Psi)) nlags = 1
Psi.new=array(0,c(dim(Psi)[1],dim(Psi)[1],nlags+1))
Psi.new[,,1]=diag(d)
if(nlags==1) {
Psi.new[,,2]=Psi
} else
for(i in 1:nlags){ Psi.new[,,i+1]=Psi[,,i]}

Psi.new=timedom(Psi.new,lags=0:nlags)

innovs =  mvtnorm::rmvnorm(n = n+nlags,mean = rep(0,d), sigma = sigma)
if (noise == 'mt'){
innovs = mvtnorm::rmvt(n = n+nlags, sigma = sigma, df=df)
}

out=filter.process(innovs,Psi.new)
out=out[(nlags+1):(n+nlags),]
out
}


## Try the freqdom package in your browser

Any scripts or data that you put into this service are public.

freqdom documentation built on Oct. 4, 2022, 5:05 p.m.