R/neighbourmatrix.R

In Sondre91/STGARCH: STGARCH

#' Create neighbourmatrix
#'
#' Function that creates a standard neighbourmatrix for spatial data observed on a equidistant grid with the circular assumption.
#' @param d Number of rows. If \code{d=1} the function creates a 1D neighbourmatrix.
#' @param M Total number of data.
#' @return Neighbourmatrix of dimension \code{M x M}.
#' @export
neighbourmatrix<-function(d,M){
  #d is one of the dimensions, M is total number of data
  d1<-d
  d2<-M/d1
  if(d2%%1!=0) stop("Dimension not a multiple of M")
  if(d1==1){
    W<-diag(M)+diag(M+1)[-c(1),-c(M+1)]+diag(M+1)[-c(M+1),-c(1)]
    W[M,1]<-W[1,M]<-1 #Adding circular corners
  }else{
    #Creating blocks
    w<-numeric(length=d1^2)
    w[1:2]<-1
    w[(d1^2-1):d1^2]<-1
    w[c(d1,d1*(d1-1)+1)]<-1
    for(i in seq(d1+1,d1*(d1-1)-2,d1+1))
      w[i:(i+2)]<-1
    W1<-matrix(w, ncol=d1,nrow=d1)
    #Placing blocks
    W<-matrix(0,ncol=prod(d1,d2),nrow=prod(d1,d2))
    for(i in 1:d2){
      W[1:d1+d1*(i-1),1:d1+d1*(i-1)]<-W1
    }
    for(i in 1:(d2-1)){
      W[1:d1+d1*i,1:d1+d1*(i-1)]<-W1
      W[1:d1+d1*(i-1),1:d1+d1*i]<-W1
    }
    W[prod(d1,d2-1)+1:d1,1:d1]<-W[1:d1,prod(d1,d2-1)+1:d1]<-W1
    #apply(W,1,sum)
    #apply(W,2,sum)

  }
  return(W)
}

#' Create reflecting neighbourmatrix
#'
#' Function that creates a reflecting neighbourmatrix for spatial data observed on a equidistant grid.
#' @param d Number of rows. If \code{d=1} the function creates a 1D neighbourmatrix.
#' @param M Total number of data.
#' @return Reflecting neighbourmatrix of dimension \code{M x M}.
#' @export
reflecting_neighbourmatrix<-function(d,M){
  #d is one of the dimensions, M is total number of data
  d1<-d
  d2<-M/d1
  if(d2%%1!=0) stop("Dimension not a multiple of M")
  if(d1==1){
    W<-diag(M)+diag(M+1)[-c(1),-c(M+1)]+diag(M+1)[-c(M+1),-c(1)]
    W[1,2]<-W[M,M-1]<-2 #Adding reflecting corners
  }else{
    #Creating blocks
    w<-numeric(length=d1^2)
    w[1:2]<-1
    w[(d1^2-1):d1^2]<-1
    #w[c(d1,d1*(d1-1)+1)]<-1
    for(i in seq(d1+1,d1*(d1-1)-2,d1+1))
      w[i:(i+2)]<-1
    w[d1+1]<-w[d1^2-d1]<-2
    W1<-matrix(w, ncol=d1,nrow=d1)
    #Placing blocks
    W<-matrix(0,ncol=prod(d1,d2),nrow=prod(d1,d2))
    for(i in 1:d2){
      W[1:d1+d1*(i-1),1:d1+d1*(i-1)]<-W1
    }
    for(i in 1:(d2-1)){
      W[1:d1+d1*i,1:d1+d1*(i-1)]<-W1
      W[1:d1+d1*(i-1),1:d1+d1*i]<-W1
    }
    W[1:d1,d1+1:d1]<-W[d1^2-d1+1:d1,d1^2-2*d1+1:d1]<-2*W[1:d1,d1+1:d1]
    #apply(W,1,sum)
    #apply(W,2,sum)

  }
  return(W)
}