以下是部分代码

rm(list=ls())

library(tseries)   
library(fGarch)
#library(rugarch)
library(zoo)

source("ind_test.R")
source("bern_test.R")

# data
p = get.hist.quote(instrument = "^djc", start = "2000-01-30",end="2017-11-01",quote="Close",quiet=T)
y=diff(log(p))      
y=coredata(y)


# Set backtest up
TT = length(y)      
WE = 500          # est. window
p = 0.05          # probability level
l1 = WE*p         # location of "-VaR(p)" in a WE-dimensional distribution
value = 1;        # portfolio value

# initialize forecasts matrix (will use 4 models)
VaR = matrix(nrow=TT,ncol=4)  

# EWMA parameters
lambda = 0.94;
s11 = var(y[1:30]);    # initialisation of sigma^2
for(t in 2:WE) s11 = lambda * s11 + (1-lambda) * y[t-1]^2 # this stops at WE

# Running the backtest
for (t in (WE+1):TT){
  t1 = t-WE;         # start of the data window
  t2 = t-1;          # end of the data window
  window = y[t1:t2]  # data for estimation
  # ewma
  s11   = lambda * s11  + (1-lambda) * y[t-1]^2 # this starts at WE+1
  VaR[t,1] = -qnorm(p) * sqrt(s11) * value
  # ma
  VaR[t,2] = - qnorm(p)* sd(window) * value
  # hs
  ys = sort(window)         
  VaR[t,3] = -ys[l1]*value  
  
  # garch(1,1)
  g=garchFit(formula = ~ garch (1,1), window ,trace=FALSE, include.mean=FALSE)
  par=g@fit$matcoef        
  s4=par[1]+par[2]*window[WE]^2+par[3]*g@h.t[WE]   # h.t is cond variance
  VaR[t,4] = -qnorm(p) * sqrt(s4) * value # garch(1,1) VaR
}