Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(pracma)
- library(cubature)
- fYth=function(x,c1,sig) {
- y=x[1]
- th=x[2]
- res= (
- 2*exp((c1^2*(1 - th + 2*th^2))/(2*sig^2*(-1 + th)) + (y*((2*c1)/(-1 + th) + y/th))/(2*sig^2) -
- (c1^2*th^3 - c1*(-2 + th)*th*y + (-1 + th)*y^2)/(sig^2*(-1 + th)*th))*sig*y-
- (c1*exp((c1^2*th*(1 + th))/(2*sig^2*(-1 + th)) + (y*((2*c1)/(-1 + th) + y/th))/(2*sig^2) - (c1^2*th^3 - c1*(-2 + th)*th*y + (-1 + th)*y^2)/
- (sig^2*(-1 + th)*th))*sqrt(2*pi)*sqrt(1 - th)*y)+
- c1*exp(-((c1^2*th^3 - c1*(-2 + th)*th*y + (-1 + th)*y^2)/(sig^2*(-1 + th)*th)) + ((2*c1*th^2*(c1 + y))/(-1 + th) + (-(c1*th) + y)^2)/
- (2*sig^2*th))*sqrt(2*pi)*sqrt(1 - th)*y*erf((c1*sqrt(1 - th))/(sqrt(2)*sig))
- )/(2*pi*sig^3*th*sqrt((1 - th)*th))
- return(res)
- }
- FY=function(y,c1,sig) {
- a1=adaptIntegrate(fYth,lowerLimit=c(0,0.000001),upperLimit =c(y,0.999999),
- c1=c1,sig=sig,maxEval=10000)
- if (is.na(a1$integral) | a1$integral>1) print(c(y,a1$returnCode,a1$integral))
- return(a1$integral)
- }
- n=1000
- nsim=10000
- mu=0.5
- sig=2
- x=matrix(rnorm(n*nsim,mu,sig),ncol=n)
- y=apply(x,1,function(x1) max(cumsum(x1)))
- plot(ecdf(y))
- pts=FYs=c(350:700)
- for (i in 1:length(pts)) FYs[i]=FY(pts[i],mu*n,sig*sqrt(n))
- lines(pts[!is.na(FYs)],FYs[!is.na(FYs)],col="red")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement