1

我正在尝试在 R 中以数值方式计算双积分的值。我的代码和我尝试编写的方程式如下：

Expectedrsquare<-function(theta=.05,sigma=.0135,alpha=0.05,r0=0,Maturity=20){

library(mvtnorm)

DoubleIntegration<-integrate(function(s) {

sapply(s, function(s) {

  integrate(function(u){

    Mus<- exp(-alpha*s)*r0+(1-exp(-alpha*s))*theta

    Sigmas<-((sigma^2/(2*alpha))*(1-exp(-2*alpha*s)))^0.5

    Muu<- exp(-alpha*u)*r0+(1-exp(-alpha*u))*theta

    Sigmau<-((sigma^2/(2*alpha))*(1-exp(-2*alpha*u)))^0.5

    Cov<-(sigma^2/(2*alpha))*exp(-alpha*(s-u))*(1-exp(-2*alpha*u))

    Chi<-Cov/(Sigmas*Sigmau)

    Zetas<-Mus/Sigmas

    Zetau<-Muu/Sigmau
         Phi2<-pmvnorm(lower=-Inf*c(1,1),upper=c(-Zetas,-Zetau),mean=c(0,0),corr=rbind(c(1,Chi),c(Chi,1)))

    IntegralFunction<-(Mus*Muu+Cov)*(1-pnorm(-Zetas)-pnorm(-Zetau)-Phi2[1])+
      Sigmas*Muu*dnorm(Zetas)*pnorm((Zetau-Chi*Zetas)/((1-Chi^2)^0.5))+
      Sigmau*Mus*dnorm(Zetau)*pnorm((Zetas-Chi*Zetau)/((1-Chi^2)^0.5))+
      Sigmas*Sigmau*(((1-Chi^2)/(2*pi))^0.5)*dnorm(((Zetau^2-2*Chi*Zetas*Zetau+Zetas^2)/((1-Chi^2)))^0.5)

    return(IntegralFunction)

    }, 0, s)$value

  })

}, 0, Maturity)

return(2*DoubleIntegration)
}

方程

当我编译并运行代码时，出现以下错误：

Error in checkmvArgs(lower = lower, upper = upper, mean = mean, corr = corr,  : 
  ‘diag(corr)’ and ‘lower’ are of different length

有人可以帮我解决错误或建议另一种解决问题的方法。

谢谢。

4

2 回答 2

0

错误是明确的：

相关性和较低的参数应该具有相同的长度。

lower  =-Inf*c(1,1)
correl = rbind(c(1,Chi),c(Chi,1)

我认为问题来自Chi的计算。

的预期维度是Sigmau多少？Sigmas?? 我猜Sigmas是问题所在。

于 2013-11-26T23:02:10.250 回答
0

终于得到了这个工作。只需在循环中分别计算双变量 CDF 并将其矢量化。

Zetas<-Mus/Sigmas
Zetau<-Muu/Sigmau
L<-length(u)
Phi2<-vector(length=L)
for(i in 1:L){
     t<-pmvnorm(lower=-Inf*c(1,1),upper=c(-Zetas,-Zetau[i]),mean=c(0,0),
                corr=rbind(c(1,Chi[i]),c(Chi[i],1)))
      Phi2[i]<-t[1]
    }
IntegralFunction <- (Mus*Muu+Cov)*(1-pnorm(-Zetas)-pnorm(-Zetau)-Phi2)+
      Sigmas*Muu*dnorm(Zetas)*pnorm((Zetau-Chi*Zetas)/((1-Chi^2)^0.5))+
      Sigmau*Mus*dnorm(Zetau)*pnorm((Zetas-Chi*Zetau)/((1-Chi^2)^0.5))+
      Sigmas*Sigmau*(((1-Chi^2)/(2*pi))^0.5)*
                    dnorm(((Zetau^2-2*Chi*Zetas*Zetau+Zetas^2)/((1-Chi^2)))^0.5)
于 2013-12-03T16:11:00.227 回答