In clulu1014/StatComp21024: Stochastic Block Model (SBM) for network vector autoregression (NAR(1)) structure

knitr::opts_chunk$set(echo = TRUE)

Overview

StatComp21024 is a R package developed to Simulate NAR(1) model.

\

Consider a large-scale social network with $N$ nodes indexed by $1\leq i\leq N$.

Define an adjacency matrix $A=(a_{i_1i_2})\in\mathbb{R}^{N\times N}$, where $a_{i_1i_2} = 1$ if there exists a social relationship from $i_1$ to $i_2$, and $a_{i_1i_2} = 0$ otherwise. And do not allow any node to be self-related, i.e., $a_{ii}=0$ for any $1\leq i\leq N$. We assume $A^T=A$.

Let $Y_{it}\in\mathbb{R}^1$be the continuous response obtained from node $i$ at time point $t$.

\

Under a network framework, $Y_{it}$ might be affected by four different factors:

1.momentum effect: $Y_{i(t-1)}$

2.network effect: $n^{-1}i\sum_ja{ij}Y_{j(t-1)}$

3.nodal effect: $Z_i$

4.an independent noise: $\varepsilon_{it}$

\ \

The NAR(1) model: $$Y_{it}=\beta_1\sum\limits_{j=1}^Na_{ij}Y_{j(t-1)}+\beta_2Y_{i(t-1)}+\varepsilon_{it},~1\leq i\leq N,~1\leq t\leq T$$

$$\mathbb{Y}{t}=\mathbb{X}{t-1}\beta+\varepsilon_t,~\mathbb{X}t=(X{1t},...,X_{Nt}),~X_{1t}=(w_i^T\mathbb{Y}t,Y{it}),~w_{i}=(a_{ij}/n_i:1\leq j\leq N)^T,~\beta=(\beta_1,\beta_2)^T$$ \

The netA function is used to generate adjacency matrix A.

netA<- function(K,N){
  A<- matrix(NaN,N,N) 
  block<- sample(1:K,size=N,replace = T,prob=rep(1/K,K)) 
  for (i in 1:N){
    for (j in 1:i){
      if(block[i]==block[j]) 
      { A[i,j]<-A[j,i]<- sample(c(0,1),size=1,replace = T,prob=c(0.2,0.8))  } 
      else 
      { A[i,j]<-A[j,i]<- sample(c(0,1),size=1,replace = T,prob=c(0.9,0.1)) }
    }
  }
  diag(A)<-0 
  A<- A[which(rowSums(A)!=0),which(rowSums(A)!=0)]
  return(A)
}

\

Example of netA function:

We assume there are 3 blocks and the network size is 10. Use the netA function to generate a Stochastic Block Model $A$.

set.seed(1234)
library(StatComp21024)
K<- 3
N<- 10
A<- netA(K,N)
A

The obR function is used to generate time series $Y_t$ by the adjacency matrix $A$.

obR<- function(A,beta,T){ 
  N<- nrow(A)
  e<- matrix(rnorm(N*T,0,1),N,T)
  W<- diag(1/rowSums(A)) %*% A
  G<- beta[1]*W + beta[2]*diag(N)
  Y<- matrix(0,N,T)
  X<- list()
  for (s in 1:T){
    X[[s]]<- matrix(0,N,2)
  }
  Y[,1]<- rnorm(N,0,1)
  for (i in 2:T){
    Y[,i]<- G %*% Y[,i-1]+e[,i]
  }
  for (j in 1:T){
    for (k in 1:N){
      X[[j]][k,]<- rbind(W[k,] %*% Y[,j], Y[k,j])
    }
  }
  return(list(Y=Y,X=X))
}

\ Example of obR function:

According to the NAR(1) model, we assume the coefficient vector is $\beta=(0.2,-0.4)^T$, and randomly generate the sample $(\mathbb{Y}_t,\mathbb{X}_t)$.

set.seed(1234)
library(StatComp21024)
T<- 100
K<- 3
N<- 10
beta<- c(0.2,-0.4)
A<- netA(K,N)
obR(A,beta,T)

\

The betahatR function uses the $(\mathbb{Y}_t,\mathbb{X}_t)$ which generated by obR function and estimates $\beta$ by

$$\hat{\beta}=\Big(\sum\limits_{t=1}^T\mathbb{X}^T_{t-1}\mathbb{X}{t-1}\Big)^{-1}\sum\limits{t=1}^T\mathbb{X}^T_{t-1}\mathbb{Y}_t$$

betahatR<- function(Y,X){
  Z<-V<- list()
  SumZ<- matrix(0,2,2)
  SumV<- matrix(0,2,1)
  for (t in 2:ncol(Y)){
    Z[[t-1]]<- t(X[[t-1]]) %*% X[[t-1]]
    SumZ<- Z[[t-1]]+SumZ
    V[[t-1]]<- t(X[[t-1]]) %*% Y[,t]
    SumV<- V[[t-1]]+SumV
  }
  betahat<- solve(SumZ) %*% SumV
  return(list(betahat=betahat,ni=solve(SumZ)))
}

\

Example of betahatR function:

Use the OLS to estimate the coefficients.

set.seed(1234)
library(StatComp21024)
T<- c(10,50,100,1e4)
K<- 3
N<- 10
beta<- c(0.2,-0.4)
A<- netA(K,N)
res<- list()
hatbeta<- matrix(0,2,length(T))
for (t in 1:length(T)){
res[[t]]<- obR(A,beta,T[t])
hatbeta[,t]<- betahatR(res[[t]]$Y,res[[t]]$X)$betahat
}
rownames(hatbeta)<- c('beta1','beta2')
colnames(hatbeta)<- T
hatbeta

The estimated value gets closer to the true value as the sample size ($T$) increases.

\

The RMSE function simulates the following measures to gauge estimators' performances.

The root mean square error is evaluated by $$RMSE_j={R^{-1}\sum_{r=1}^R(\hat{\beta}_j-\beta_j)^2}^{1/2}$$

Next, for each replication, a 95% confidence interval is constructed for $\beta_j$ as $$CI_j^{(r)}=(\hat{\beta}j^{(r)}-z{0.975}\widehat{SE}j^{(r)},\hat{\beta}_j^{(r)}+z{0.975}\widehat{SE}j^{(r)})$$ where $\widehat{SE}_j^{(r)}$ is root square of the $j$ th diagonal element of $\Big(\sum\limits{t=1}^T\mathbb{X}^T_{t-1}\mathbb{X}{t-1}\Big)^{-1}\hat{\sigma}^2$ with $\hat{\sigma}^2=(NT)^{-1}\sum{i,t}(Y_{it}-X_{it}^T\hat{\beta}^{(r)})^2$. $z_{\alpha}$ is the $\alpha$th quantile of a standard normal distribution.

The coverage probability (CP) is computed as

$$CP_j=R^{-1}\sum_{r=1}^RI(\beta_j\in CI_j^{(r)})$$ where $I(\cdot)$ is the indicator function.

RMSE<- function(rep,K,N,beta,T){
  hatbeta<-sigmahat<- CIl<- CIr<- matrix(0,2,rep)
  eps<- numeric(T)
  for (r in 1:rep){
    A<- netA(K,N)
    res<- obR(A,beta,T)
    hatbeta[,r]<- betahatR(res$Y,res$X)$betahat
    for (t in 1:T){
      eps[t]<- t(res$Y[t]-res$X[[t]] %*% hatbeta[,r]) %*% (res$Y[t]-res$X[[t]] %*% hatbeta[,r])
    }
    sigmahat[,r]<- sqrt(diag(betahatR(res$Y,res$X)$ni*sum(eps)/(T*nrow(A))))
    CIl[,r]<- hatbeta[,r]-qnorm(0.975)*sigmahat[,r]
    CIr[,r]<- hatbeta[,r]+qnorm(0.975)*sigmahat[,r]
  }
  RMSE1<- (sum((hatbeta[1,]-beta[1])^2)/rep)^(1/2)
  RMSE2<- (sum((hatbeta[2,]-beta[2])^2)/rep)^(1/2)
  CPbeta1<- mean(CIl[1,]<=beta[1] & beta[1]<=CIr[1,])
  CPbeta2<- mean(CIl[2,]<=beta[2] & beta[2]<=CIr[2,])
  return(rbind(CPbeta1,CPbeta2,RMSE1,RMSE2))
}

\ Example of RMSE function:

set.seed(1234)
library(StatComp21024)
rep<- 1000
T<- 100
K<- 3
N<- 10
beta<- c(0.2,-0.4)
RMSE(rep,K,N,beta,T)

clulu1014/StatComp21024 documentation built on Dec. 23, 2021, 11:14 p.m.