在 R 中结合 3 个不同的矩阵图

Question

谁能告诉我如何创建一个包含 3 个不同矩阵数据集的图。一般来说，我有 3 个不同的数据矩阵，所有 1*1001 维度，我希望将所有 3 个数据都绘制在同一张图上。

我设法一次绘制一个矩阵，并组装代码以创建其他 2 个矩阵，但不绘制它。 B[i,] 是随机生成的数据。我想知道的是，将所有 3 个图放在一张图上的编码是什么。

一个矩阵的代码： ntime

for(i in 1:(ntime+1)){
average.price.at.each.timestep[i]<-mean(B[i,])
}

matplot(t, t(average.price.at.each.timestep), type="l", lty=1, main="MC Price of a Zero Coupon Bond", ylab="Price", xlab = "Option Exercise Date") 

3 的代码：

average.price.at.each.timestep<-matrix(0,nrow=1,ncol=ntime+1)
s.e.at.each.time <-matrix(0,nrow=1,ncol=ntime+1)
upper.c.l.at <- matrix(0,nrow=1,ncol=ntime+1)
lower.c.l.at <- matrix(0,nrow=1,ncol=ntime+1)
std <- function(x) sd(x)/sqrt(length(x))

for(i in 1:(ntime+1)){
average.price.at.each.timestep[i]<-mean(B[i,])
s.e.at.each.time[i] <- std(B[i,])
upper.c.l.at[i] <- average.price.at.each.timestep[i]+1.96*s.e.at.each.time[i]
lower.c.l.at[i] <- average.price.at.each.timestep[i]-1.96*s.e.at.each.time[i]
}

我仍在为此苦苦挣扎，因为我无法获得与我的数据集相匹配的解决方案，我现在将生成矩阵 B 的以下代码作为一个工作示例包含在内，这样您就可以看到我正在处理的数据.正如您所看到的，它生成了一个不同价格的图，我想要一个包含平均价格和平均值置信区间的图。

# Define Bond Price Parameters
#
P<-1                            #par value

# Define Vasicek Model Parameters
#
rev.rate<-0.3                   #speed of reversion
long.term.mean<-0.1             #long term level of the mean
sigma<-0.05                  #volatility
r0<-0.03                            #spot interest rate
Strike<-0.05
# Define Simulation Parameters
#
T<-50                           #time to expiry
ntime<-1000                 #number of timesteps
yearstep<-ntime/T               #yearstep
npaths<-1000                    #number of paths
dt<-T/ntime                     #timestep
R <- matrix(0,nrow=ntime+1,ncol=npaths) #matrix of Vasicek interest rate values
B <- matrix(0,nrow=ntime+1,ncol=npaths) # matrix of Bond Prices 

R[1,]<-r0                   #specifies that all paths start at specified spot rate
B[1,]<-P

# do loop which generates values to fill matrix R with multiple paths of Interest Rates as they evolve over time.
# stochastic process based on standard normal distribution

for (j in 1:npaths) {
 for (i in 1:ntime) {
   dZ <-rnorm(1,mean=0,sd=1)*sqrt(dt)
   Rij<-R[i,j]
   Bij<-B[i,j]
   dr <-rev.rate*(long.term.mean-Rij)*dt+sigma*dZ
   R[i+1,j]<-Rij+dr
   B[i+1,j]<-Bij*exp(-R[i+1,j]*dt)   
  }
} 

t<-seq(0,T,dt)
par(mfcol = c(3,3))


matplot(t, B[,1:pmin(20,npaths)], type="l", lty=1, main="Price of a Zero Coupon Bond", ylab="Price", xlab = "Time to Expiry") 

Answer 1

您的示例无法重现，因此我创建了一些假数据，希望其结构与您的相似。如果这不是您想要的，请告诉我，我会根据需要进行更新。

# Fake data
ntime <- 100 
mat1 <- matrix(rnorm(ntime+1, 10, 2), nrow=1, ncol=ntime+1)
mat2 <- matrix(rnorm(ntime+1, 20, 2), nrow=1, ncol=ntime+1)
mat3 <- matrix(rnorm(ntime+1, 30, 2), nrow=1, ncol=ntime+1)

matplot(1:(ntime+1), t(mat1), type="l", lty=1, ylim=c(0, max(c(mat1,mat2,mat3))),
        main="MC Price of a Zero Coupon Bond", 
        ylab="Price", xlab = "Option Exercise Date") 

# Add lines for mat2 and mat3
lines(1:101, mat2, col="red")
lines(1:101, mat3, col="blue")

更新：这是你想要做的吗？

matplot(t, t(average.price.at.each.timestep), type="l", lty=1, 
        main="MC Price of a Zero Coupon Bond", ylab="Price", 
        xlab = "Option Exercise Date") 
matlines(t, t(upper.c.l.at), lty=2, col="red")
matlines(t, t(lower.c.l.at), lty=2, col="green")

见下图。如果您有多个要绘制的列（如在更新的示例中绘制 20 个单独的路径）并且您想为所有列添加较低和较高的 CI（尽管这会使绘图不可读），只需使用矩阵与average.price.at.each.timestep 中的每个路径对应的CI 值的上限和下限，并使用matlines 将它们添加到现有的多路径图中。

Answer 2

使用ggplot2 和reshape2 可以做到这一点。您拥有的结构有点笨拙，您可以通过使用数据框而不是矩阵来改进。

#Dummy data
average.price.at.each.timestep <- rnorm(1000, sd=0.01)
s.e.at.each.time <- rnorm(1000, sd=0.0005, mean=1)

#CIs (note you can vectorise this):
upper.c.l.at <- average.price.at.each.timestep+1.96*s.e.at.each.time
lower.c.l.at <- average.price.at.each.timestep-1.96*s.e.at.each.time

#create a data frame:
prices <- data.frame(time = 1:length(average.price.at.each.timestep), price=average.price.at.each.timestep, upperCI= upper.c.l.at, lowerCI= lower.c.l.at)

library(reshape2)
#turn the data frame into time, variable, value triplets
prices.t <- melt(prices, id.vars=c("time"))

#plot
library(ggplot2)
ggplot(prices.t, aes(time, value, colour=variable)) + geom_line()

这会产生以下情节：

这可以通过使用geom_ribbon 来改善：

ggplot(prices, aes(time, price)) + geom_ribbon(aes(ymin=lowerCI, ymax=upperCI), alpha=0.1) + geom_line()

这会产生这个情节：

Answer 3

这是另一个稍有不同的 ggplot 解决方案，它不需要您先计算置信限 - ggplot 会为您完成。

# create sample dataset
set.seed(1)   # for reproducible example
B <- matrix(rnorm(1000,mean=rep(10+1:10/2,each=10)),nc=10)

library(ggplot2)
library(reshape2)   # for melt(...)
gg <- melt(data.frame(date=1:nrow(B),B), id="date")
ggplot(gg, aes(x=date,y=value)) + 
  stat_summary(fun.y = mean, geom="line")+
  stat_summary(fun.y = function(y)mean(y)-1.96*sd(y)/sqrt(length(y)), geom="line",linetype="dotted", color="blue")+
  stat_summary(fun.y = function(y)mean(y)+1.96*sd(y)/sqrt(length(y)), geom="line",linetype="dotted", color="blue")+
  theme_bw()

stat_summary(...) 总结了给定 x 值（日期）的 y 值。所以在第一次调用中，它计算平均值，第二次计算 lowerCL，第三次计算 upperCL。

您还可以创建一个CL(...) 函数，并调用它：

CL <- function(x,level=0.95,type=c("lower","upper")) {
  fact <- c(lower=-1,upper=1)
  mean(x) - fact[type]*qnorm((1-level)/2)*sd(x)/sqrt(length(x))
}
ggplot(gg, aes(x=date,y=value)) + 
  stat_summary(fun.y = mean, geom="line")+
  stat_summary(fun.y = CL, type="lower", geom="line",linetype="dotted", color="blue")+
  stat_summary(fun.y = CL, type="upper", geom="line",linetype="dotted", color="blue")+
  theme_bw()

这会产生一个与上面相同的图。