金融计量经济学中使用R构建置信区间

最近有很多同学问我如何构建置信区间，所以简单介绍一下学习成果

一般来说要我们构建置信区间的呢，是未知的数，我们得到的都是计算值估计值或者预测值，比如说计算一个集合X的均值的置信区间，计算模型参数的置信区间，以及计算预测值的置信区间。

那么如何计算置信区间的，我们可以计算很多次该值，就可以得到它的分布了，根据大数定律，我们得到的分布就是真实值的分布。

但是如果我们使用的数据集是同一个数据集，那我们计算出的值是不会有变化的，那么如何让他有变化，但是又保证得到的结果是可信的结果呢。

一个方法就是蒙特卡洛模拟，比如说我们有一个集合X，那么我们可以计算它的均值和方差，注意这里我们需要假设集合X是服从正态分布的，然后我们随机生成一个和X一样分布的集合，用这样的方法我们可以生成N个和X一样分布的不同的集合，那么我们利用这些集合计算出N个不同的均值，这些均值放在一起，就是真实值的分布。N无限大的时候，这个均值集合的分布就是正态分布，当然我们做不到无穷大，但我们把它当成正态分布，所以我们在寻找置信区间的上下限时，只需要寻找到百分位数，例如我们模拟了2000次，90%置信区间的上限就是将我们得到的2000个均值从小到大排的第1900个，下限就是第100个这样子。（请不要都举这个例子）

#下面是均值的例子
#这里假设我的X里已经存了数据了
#样本量为：
N <- NROW(X)
# X的均值为
Xmean <- mean(X)
# X的方差为
Xsd <- sd(X)
K <- 1000 #模拟次数
result <- matrix(NA, R, 1) #创建一个空矩阵来存模拟结果
#然后我们模拟K个X
for(i in 1:K){ #循环做这件事K次
	xsim <- rnorm(N, mean=Xmean, sd=Xsd) #模拟一个和X同分布的大小为N的样本
	result[i,1] <- mean(xsim) #将这个样本的均值存在result里面
}
#求95%的置信区间
quantile(result, probs=c(2.5/100, 97.5/100)

那么如果是模型呢，例如y=a+bx+error，这种情况下，我们首先是得对数据集跑出a和b，然后就像上面的例子一样，模拟X，再根据y=a+bx+error模拟出Y，同样做N个模拟，跑出N个不同的a和b，再视你对哪个感兴趣就对哪个做置信区间。那么如果是预测，就先模拟出数据，然后拟合模型，然后得到预测值，这样就可以得到很多个不同的预测值。

#下面是AR（1）模型的例子,假如我们已经有时间序列X
#首先我们要拟合一个AR模型y_t=a+b*y_(t-1)+error
est <- arima(X, order=c(1,0,0))
arcoefs <- coeftest(est) #取出参数的测量值
a <- arcoefs[2,1]
b <- arcoefs[2,2]
K <- 1000 #模拟1000次，你也可以模拟更多或更少次，但不要太少
result <- matrix(NA, R, 2) #创建一个空矩阵来存a和b的模拟值
colnames <- c("b", "a")#第一列存b第二列存a
for(i in 1:K){
	xsim <- arima.sim(list(order=c(1,0,0), ar=b), n=N) #模拟b*y_(t-1)+error
	xsim <- a+xsim #模拟y_t
	
	#用xsim来拟合AR（1）
	est.sim <- arima(xsim, order=c(1,0,0))
	
	#把结果存起来
	coef.sim <- est.sim$coef
	result[i,] <- coef.sim
}

#计算95%的置信区间
#a的置信区间
quantile(result[,2], probs=c(2.5/100, 97.5/100)
#b的置信区间
quantile(result[,1], probs=c(2.5/100, 97.5/100)

其实bootstrap和Monte Carlo的原理是很像的，只是蒙特卡洛的数据是通过模拟出来的，因此需要假设数据服从特定分布，但是bootstrap使用了一种聪明的方法，它通过抽样来模拟样本，这样既做到了不改变样本的分布，又不需要对样本的分布进行假设。
那么主要有IID bootstrap，就是有放回的抽取样本组成新样本；stock bootstrap，主要用于时间序列，就是把样本装到block里面，每个block里面的数据是有顺序的，然后从block里面抽取block组成新的样本；还有residual bootstrap，用于在存在单位根的时候，但是一般来说我们在对数据进行模型拟合之前，第一步是让数据平稳，所以对最后一种就不多做介绍了。

#首先说IID bootstrap，假设我们已经有一个平稳的数据X
#样本大小
N <- NROW(X)
#构建随机因子,1-N之间的N个随机数
select <- round(runif(n=N, min=1, max=N)
x_new <- X[select]

这就是一次iid bootstrap，将该过程像蒙特卡洛一样进行N次得到N个新的样本，并以此对要构建置信区间的值进行N次计算，可以得到该值的近似分布，因为过程大同小异就不写了（懒）

#下面是block bootstrap，假设我们已经有一个平稳的数据X
#样本大小
N <- NROW(X)
#假设block的大小是b，是一个已赋值的变量
x_new <- NULL
k <- NROW(x_new) #就是0
while(k<=(N+5){#保证新样本的样本数不比本来的样本少
	start <- runif(1,1,N-b+1) #生成以1到N-b+1之间的随机数，作为block的起始数据定位
	end <- start+b-1 #block的终止数据定位
	block <- X[(start+1):end] #选择出的block
	x_new <- c(x_new, block) #将选出的block加进去
	k <- NROW(x_new) # k <- k+b
}
x_new <- x_new[1:N] #考虑到这样可能新的样本数量会多于原样本，只保留前N个数据