（随机微分方程）SDE模拟股价变化

library("readr")
library("dplyr")
library(PerformanceAnalytics)
setwd("e:/R/tail risk/day/sp")
# 利用GBM模型模拟万科A股价变化
sp1<-read.csv("sp1.csv")%>%select(1,2,5)
wka<-filter(sp1,Stkcd==2&Trdmnt<=20181231&Trdmnt>=20180101)%>%select(-1)%>%arrange(Trdmnt)
wka1<-mutate(wka,ln1=log(Clsprc))
wka2<-bind_cols(slice(wka1,-n()),select(slice(wka1,-1),3))%>%mutate(ln=ln11-ln1)
#计算股价样本均值及标准差
a<-mean(wka2$ln);b<-sd(wka2$ln)
#估计出mu与sigma的值
mu<-a+b^2/2;sigma<-b# 可得出mu=-0.0009315989，sigma=0.02683159
##代入GBM 模型进行模拟可得出股价
library(gcookbook)
library(ggplot2)
library(reshape2)
for(i in 2:nrow(wka)){
  wka[i,3]=mean(wka[i-1,2]*exp(sigma*(rnorm(1000,0,1))+(mu-sigma^2/2)))
}
wka<-rename(wka,real=Clsprc,simulate=V3)%>%na.omit()
wka3<-melt(wka,id.vars="Trdmnt")%>%rename(price=value)
ggplot(wka3,aes(x=factor(Trdmnt),y=price,group=1))+geom_line(aes(color=variable))+labs(x="Time")
plot(wka[,3]-wka[,2],type="l",main="残差图",ylab="pricespread",xlab="day")
#可知拟合优度为：
Q<-1-sum(((wka[,3]-wka[,2])^2)[-1])/sum(wka[,2]^2) # 拟合优度Q=0.99927

注：红色是2018年万科A的真实股价、绿色是利用布朗运动模型估计的股价。