R/hofa.sim.R

In GuanglinHuang/hofa: higher-order multi-cumulant factor analysis

#' Generate panel data from [H]igher-[O]rder multi-cumulant [F]actor [A]nalysis of DGP1
#'
#' @param n An integer, the number of variables.
#' @param t An integer, the number of observations.
#' @param k An integer, the number of factors.
#' @param par_f A list of four elements, the parameters of factor distribution, the first list is sig_f = (sig_f1,sig_f2,...), the second list is eta_f = (eta_f1,eta_f2,...), the third list is p_f = (p_f1,p_f2,...), the fourth list is q_f = (q_f1,q_f2,...).
#' @param par_e A list of four elements, the parameters of error distribution, par_e = list(sig_e,eta_e,p_e,q_e).
#' @param alpha The eigenvalue rate of error terms.
#' @param rho_f A k x 1 vector, the auto-regressive parameters of factors.
#' @param rho_e The auto-regressive parameters of errors.
#' @param ... Any other parameters.
#' @return A T x N matrix, generated by DGP1 of hofa.
#' @examples
#' n = 100
#' t = 200
#' k = 2
#' par_f = list(rep(1,k),rep(0.8,k),rep(1,k),rep(Inf,k))
#' par_e = list(1,0,2,Inf)
#' rho_f = c(0.5,0.2)
#' alpha = 0
#' hofa.DGP1(n,t,k,par_f,par_e,alpha,rho_f,rho_e = 0.2)



hofa.DGP1 = function(n,t,k,par_f,par_e,alpha,rho_f,rho = 0.5,theta = 1,rho_e = 0.2,...){
  sig_f <- par_f[[1]]
  lam_f <- par_f[[2]]
  p_f <- par_f[[3]]
  q_f <- par_f[[4]]

  sig_e <- par_e[[1]]
  lam_e <- par_e[[2]]
  p_e <- par_e[[3]]
  q_e <- par_e[[4]]

  FF = matrix(NA,t,k)
  for (ii in 1:k) {
    FF[,ii] <- arima.sim(list(ar = rho_f[ii]),t,rand.gen = function(n,...){rsgt(n,0,sig_f[ii],lam_f[ii],p_f[ii],q_f[ii])})
  }

  if(rho_e != 0){
    et = matrix(0,t,n)
    for (j in 1:n) {
      et[,j] <- arima.sim(list(ar = rho_e),t,rand.gen = function(n,...){rsgt(n,0,sig_e,lam_e,p_e,q_e)})
    }
  }

  if(rho_e == 0){
    et = matrix(NA,t,n)
    for (j in 1:n){
      et[,j] = rsgt(t,0,sig_e,lam_e,p_e,q_e,var.adj = T)
    }
  }

  W  = matrix(rnorm(k*n),n,k)


  d = diag(theta*(1:n)^{-rho})

  u = pracma::randortho(n,"orthonormal")

  O = u%*%d%*%t(u)
  O_half = sqrt(d)%*%t(u)

  UU = et

  E = UU%*%O_half

  X = FF%*%t(W) + E

  return(list(X = X,FF = FF,W = W,E = E))
}


#' Generate panel data from [H]igher-[O]rder multi-cumulant [F]actor [A]nalysis of DGP2
#'
#' @param n An integer, the number of variables.
#' @param t An integer, the number of observations.
#' @param k An integer, the number of factors.
#' @param par_f A list of four elements, the parameters of factor distribution, the first list is sig_f = (sig_f1,sig_f2,...), the second list is eta_f = (eta_f1,eta_f2,...), the third list is p_f = (p_f1,p_f2,...), the fourth list is q_f = (q_f1,q_f2,...).
#' @param par_e A list of four elements, the parameters of error distribution, par_e = list(sig_e,eta_e,p_e,q_e).
#' @param par_cove A list of four elements, the parameters control the covariance structure of the error terms, par_cove = list(beta,J,rho,msig_e).
#' @param rho_f A k x 1 vector, the auto-regressive parameters of factors.
#' @param ... Any other parameters.
#' @return A T x N matrix, generated by DGP2 of hofa.
#' @examples
#' n = 100
#' t = 200
#' k = 2
#' par_f = list(rep(1,k),rep(0.8,k),rep(1,k),rep(Inf,k))
#' par_e = list(1,0,2,Inf)
#' rho_f = c(0.5,0.2)
#' par_cove = list(beta = 0.2,J = n/10,rho = 0.2,msig_e = c(1,5))
#' hofa.DGP2(n,t,k,par_f,par_e,par_cove,rho_f)

hofa.DGP2 = function(n,t,k,par_f,par_e,par_cove,rho_f,...){

  sig_f <- par_f[[1]]
  lam_f <- par_f[[2]]
  p_f <- par_f[[3]]
  q_f <- par_f[[4]]

  sig_e <- par_e[[1]]
  lam_e <- par_e[[2]]
  p_e <- par_e[[3]]
  q_e <- par_e[[4]]

  beta = par_cove[[1]]
  J = par_cove[[2]]
  rho = par_cove[[3]]
  msig_e = par_cove[[4]]


  m = min(t,n)

  W=matrix(rnorm(n*k),n,k)

  theta_e = runif(n,msig_e[1],msig_e[2])

  FF = matrix(NA,t,k)
  for (ii in 1:k) {
    FF[,ii] <- arima.sim(list(ar = rho_f[ii]),t,rand.gen = function(n,...){rsgt(n,0,sig_f[ii],lam_f[ii],p_f[ii],q_f[ii])})
  }


  if(rho != 0){
    p = n

    ut = matrix(NA,t,n)
    for (j in 1:n) {
      ut[,j] = rsgt(t,0,sig_e,lam_e,p_e,q_e,var.adj = T)
    }

    aast = c(1:p)
    et = matrix(0,t,p)
    for (j in 1:p) {
      ast = aast[j]
      et[,ast] = arima.sim(list(order = c(1,0,0), ar = rho),n = t ,rand.gen = hofa::rccsgt , ast = ast , p = p,beta = beta ,ut =ut, J = J)
    }
    et = et[,-((n+1):(n+2))]
    theta = rep(1,t)%*%t(theta_e)
    E<-sqrt(theta)*et*sqrt((1-rho^2)/(1+2*J*beta^2))
  }

  if(rho == 0){
    p = n
    ut = matrix(NA,t,n)
    for (j in 1:n) {
      ut[,j] = rsgt(t,0,sig_e,lam_e,p_e,q_e,var.adj = T)
    }

    aast = c(1:p)
    et = matrix(0,t,p)
    for (j in 1:p) {
      ast = aast[j]
      et[,ast] = hofa::rccsgt(ast = ast,p = p,beta = beta ,ut =ut, J = J)
    }
    theta = rep(1,t)%*%t(theta_e)
    E<-sqrt(theta)*et*sqrt((1-rho^2)/(1+2*J*beta^2))

  }

  X = FF%*%t(W)+E

  return(list(X = X,W = W,FF = FF,E = E))
}