Approximate asian geometric option with Heston

Question

I am trying to implement Theorem 1 from this Journal in RStudio.

The journal says the it is possible to find a approximate price of a geometric asian option in a Heston setup this way: $$X_{1cGAO}=e^{0.5\hat{\sigma}-rT+a}N(d_+)-Ke^{-rT}N(d_-)$$ where $$a=log(S_0)+0.5rT- \frac{1}{2T}\int_0^T[\int_0^t\sigma^2(s)ds]dt$$ $$\hat{\sigma}^2=\frac{1}{n^2}\sum_{j=1}^n[2(n-j)+1]\int_0^{j\frac{T}{n}}\sigma^2(s)ds$$ where $n$ goes to $\infty$ and $d_-=\frac{a-log(K)}{\hat{\sigma}}$, $d_+=d_-+\hat{\sigma}$

The trick is, according to the journal to substitute the stochastic volatility from the Heston process with its deterministic expected value: $$E[\sigma^2_t]=v_0e^{-\kappa*t}+\theta(1-e^{-\kappa t})$$

And here comes my problem:

When i am trying to implement the code in RStuido (see code below) I do not get a plausibly result. I get the price $8.592864$ and when I am using Monte Carlo simulation I get about $5.430642$

Can someone see what I am doing wrong or what I am missing?

kappa=2
theta=0.04
v0=0.04
S=100
rf=0.05
K=100
T=1

MC_Heston_Std_Asian_GEO(S,K,rf,T,v0,theta,10^6,10,kappa,0,0)


#Expected value of CIR-Prices
ExpectedVarians <- function(t) {
  exp(-kappa*t)*v0+(theta)*(1-exp(-kappa*t))
}

sqrt(ExpectedVarians(1))

#My function
AsianGEO = function(S,K,rf,T,n){

  a=log(S)+0.5*rf*T-(1/(2*T))*integrate(function(y) { 
    sapply(y, function(y) {
      integrate(ExpectedVarians, 0, y)$value
        })
      }, 0, T)$value

  sum=0
  for (j in 1:n){
    sum=sum+(2*(n-j)+1)*integrate(ExpectedVarians, 0, j*(T/n))$value
  }
  sigma_hat=sqrt(sum/n^2)
  d_minus=(a-log(K))/(sigma_hat)
  d_plus=d_minus+(sigma_hat)

  (exp(0.5*sigma_hat-rf*T+a)*pnorm(d_plus)-K*exp(-rf*T)*pnorm(d_minus))
}

AsianGEO(S,K,rf,T,1000)

Added:

MC_Heston_Std_Asian_GEO <- function(S,K,rf,T,V0,theta,n,b,kappa,sigma_vol,rho) {

  values <- c(4)
  # Start the clock!
  ptm <- proc.time()

  S_pris=0
  Price_sum_i_2=0
  Price_sum=0
  delta=(T)/(b)

  Z1=matrix(rnorm(n*b),n, b)
  Z2=matrix(rnorm(n*b),n, b)

  Z_vol=Z1
  Z=rho * Z_vol + sqrt(1-rho^2)*Z2

  for (i in 1:n) {
    x_t=log(S)
    vol_i=V0
    S_produkt = 1*S
    S_mean_GEO=0

    for (j in 1:b) {

      vol_i=V0
      x_t=x_t + (rf - 0.5*max(vol_i,0))*delta+(max(vol_i,0))**(0.5)*delta**(0.5)*Z[i,j]
      vol_i = vol_i + kappa * (theta - max(vol_i,0)) * delta + sigma_vol * (max(vol_i,0))**(0.5)*(delta)**(0.5)*Z_vol[i,j]
      S_produkt = S_produkt*exp(x_t)
    }

    S_mean_GEO=S_produkt^(1/(b+1))


    C = exp(-1*rf*T)*max(S_mean_GEO-K,0)
    Price_sum=Price_sum+C
    Price_sum_i_2=Price_sum_i_2+C^2

    if (i==n) {
      print("*STANDARD HESTON MONTE CARLO MODELLEN ASIAN GEO*")
      values[1]=Price_sum/n #PRIS
      print(paste0("Option Price: ", values[1]))
      values[2]=(Price_sum_i_2/n)-(Price_sum/n)^2 #VARIANS
      values[3]=sqrt(values[2])/sqrt(n) #STD
      print(paste0("Std on estimat: ", values[3]))
    } 
  }
  # Stop the clock
  values[4]=(proc.time() - ptm)[3]
  print(paste0("TIme: ", values[4]))

  return(list(Price=values[1],Varians=values[2],Std.dev=values[3],Time=values[4]))
}