R/estimation.R

In Sondre91/STGARCH: STGARCH

#' Estimation of 2D STARCH(1)
#'
#' Function for estimating parameters of a STARCH(1) process
#' @param X Data array of form spatial1 x spatial2 x temporal
#' @param n.iterations Integer, maximum number of iterations
#' @param return_all Returns all
#' @param W Neighbourhood matrix
#' @param var_mat Logical, determines if the covariance matrix should be computed and returned
#' @param circular Logical, determines if a circular model should be assumed
#' @param LSE Logical, determines if the Least Squares Estimate should be returned
#' @return Depending on the logical inputs. Returns a vector containing various elements.
#' @export
#' @examples
#' x <- simulation_2d_arch(param = c(.3, .1),
#'                         dim = c(5, 5, 100))
#' estim <- est(X = x,
#'              n.iterations = 30,
#'              W = neighbourmatrix(5,25),
#'              circular = TRUE)
#'
#'
est<-function(X,n.iterations=50,
              return_all=FALSE,
              W=neighbourmatrix(dim(X)[1],prod(dim(X)[1:2])),
              var_mat=FALSE,
              circular=TRUE, LSE=FALSE){
  # Regular estimation
  # Sondre H?lleland, 8.jan 2016
  # M = spatial points
  # N = temporal points
  d<-dim(X)
  M<-prod(d[1:2])
  N<-d[3]
  #set.seed(seed)
  d1<-d[1]
  d2<-M/d1
  #W<-neighbourmatrix(d1,M)
  X<-tilVektor(X)
  #Estimation step:
  mat<-array(NA_real_, dim=c(M,2,N-1))
  mat[,1,]<-1
  mat[,2,1:(N-1)]<-(W%*%X[,1:(N-1)]^2)
  #Removing regular points
  if(!circular){
    rem<-c(1:d[1],seq(d[1]+1,M-2*d[1]+1,d[1]), seq(2*d[1],M-d[1],d[1]), (M-d[1]+1):M)
    mat<-mat[-rem,,]
    X<-X[-rem,]
  }
  #LEAST SQUARES
  resultLS<-LS.estimation(X,mat, var_mat)

  #MAXIMUM LIKELIHOOD
  resultML<-max_lik(init=resultLS[1:2],
                    X,mat,
                    n.iterations=n.iterations,
                    return_all=return_all,
                    var_mat=var_mat)
  if(var_mat)
    resultML<-c(resultML$theta,resultML$var)
  if(return_all){
    return(resultML)
  }else if(LSE){
    return(c(resultML, resultLS))
  }else{return(resultML)}
}

#' Least Squares Estimation of STARCH(1)
#'
#' Function for estimating parameters of a STARCH(1) process
#' @param X Data array of form spatial1 x spatial2 x temporal
#' @param mat The design array with dimension (M,2,N-1).
#' @param var_mat Logical, determines if the covariance matrix should be computed and returned
#' @return Depending on the logical inputs. Returns a vector containing the estimates
#' (and the covariance matrix if \code{var_mat = TRUE}).
#'
LS.estimation<-function(X, mat, var_mat=FALSE){
  #X is the data, arranged as a vector
  #mat is the design array with dimension (M,2,N-1).
  #var_mat is a logic vector used for decided wether or not to return the
  #covariance matrix of the estimators.
  Y<-X^2
  d<-dim(X);N<-d[2];M<-d[1];
  #Creates a temporary array in order to use the apply function
  tmp.mat<-array(NA_real_, dim=c(M,3,N-1))
  tmp.mat[,1:2,]<-mat
  tmp.mat[,3,]<-Y[,2:N]
  v<-apply(apply(tmp.mat, 3, function(x) t(x[,1:2])%*%x[,3]),1,sum)
  design<-matrix(apply(apply(mat,3,function(x)c(t(x)%*%x)),1,sum), ncol=2, nrow=2)
  #Calculate theta
  theta<-solve(design,v)
  if(!var_mat)
    return(theta)
  else{
    #Covariance matrix calculations
    var<-2*solve(design)%*%matrix(apply(apply(mat, 3,function(x) {
      return(c(t(x)%*%diag(c((x%*%theta)^2))%*%x))}),1,sum),
      ncol=length(theta), nrow=length(theta))%*%solve(design)
    return(c(theta,var))
  }
}

#' Maximum Likelihood Estimation of STARCH(1)
#'
#' Function for estimating parameters of a STARCH(1) process
#' @param init Initial parameter vector
#' @param X Data array of form spatial1 x spatial2 x temporal
#' @param mat The design array with dimension (M,2,N-1).
#' @param n.iterations Maximum number of iterations
#' @param return_all Returns all
#' @param var_mat Logical, determines if the covariance matrix should be computed and returned
#' @return Depending on the logical inputs. Returns a vector containing various elements.
#'

max_lik<-function(init, X,mat, n.iterations=20, return_all=FALSE, var_mat=FALSE){
  # init is initial parameter vector
  # X vectorized data
  # mat design matrix
  # n.iterations is number of iterations
  # return_all decides if all iterations should be returned
  # var_mat decides if covariance matrix should be calculated and returned

  # Edit: 18. april
  # Added logical convergence stop at difference of 10^(-10)

  # Initialize theta matrix
  theta<-matrix(NA_real_, ncol=2, nrow=n.iterations+1)
  #Intial value
  theta[1,]<-init
  d<-dim(X)
  w<-array(NA_real_, dim=dim(mat))
  #Iteration process
  #for(i in 1:n.iterations){
  converge=FALSE
  zero.lim=10^(-10)
  i=1
  while(i <= n.iterations & !all(converge)){
    #Calculating weights
    w<-apply(mat,3,function(x) x%*%theta[i,])
    w<-1/(2*w^2)
    sum1<-matrix(0, ncol=2, nrow=2)
    sum2<-numeric(length=2)
    for(t in 2:d[2]){
      W<-diag(w[,t-1])
      sum1<-sum1+t(mat[,,t-1])%*%W%*%mat[,,t-1]
      sum2<-sum2+t(mat[,,t-1])%*%W%*%X[,t]^2

    }
    #Calculating next theta
    theta[i+1,]<-solve(sum1,sum2)
    converge<-abs(theta[i+1,]-theta[i,])<zero.lim
    i<-i+1
  }
  #Covariance matrix
  if(var_mat){
    retur_list<-list()
    if(return_all)
      retur_list$theta<-theta # returning all iterations
    else
      retur_list$theta<-theta[i,] # returning only last iteration
    retur_list$var<-solve(sum1) # returning covariance matrix
    return(retur_list)
  }else{
    if(return_all)
      return(theta) # returning all iterations
    else
      return(theta[i,]) # returning only last iteration
  }
}