Advertisement
cheeeesy

Untitled

Apr 24th, 2021
1,736
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 1.29 KB | None | 0 0
  1. library(pracma)
  2. library(cubature)
  3.  
  4. fYth=function(x,c1,sig) {
  5.   y=x[1]
  6.   th=x[2]
  7.   res= (
  8.     2*exp((c1^2*(1 - th + 2*th^2))/(2*sig^2*(-1 + th)) + (y*((2*c1)/(-1 + th) + y/th))/(2*sig^2) -
  9.             (c1^2*th^3 - c1*(-2 + th)*th*y + (-1 + th)*y^2)/(sig^2*(-1 + th)*th))*sig*y-
  10.       (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)/
  11.                 (sig^2*(-1 + th)*th))*sqrt(2*pi)*sqrt(1 - th)*y)+
  12.       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)/
  13.                (2*sig^2*th))*sqrt(2*pi)*sqrt(1 - th)*y*erf((c1*sqrt(1 - th))/(sqrt(2)*sig))
  14.   )/(2*pi*sig^3*th*sqrt((1 - th)*th))
  15.   return(res)
  16. }
  17.  
  18. FY=function(y,c1,sig) {
  19.   a1=adaptIntegrate(fYth,lowerLimit=c(0,0.000001),upperLimit =c(y,0.999999),
  20.                     c1=c1,sig=sig,maxEval=10000)
  21.   if (is.na(a1$integral) | a1$integral>1) print(c(y,a1$returnCode,a1$integral))
  22.   return(a1$integral)
  23. }
  24.  
  25. n=1000
  26. nsim=10000
  27. mu=0.5
  28. sig=2
  29. x=matrix(rnorm(n*nsim,mu,sig),ncol=n)
  30. y=apply(x,1,function(x1) max(cumsum(x1)))
  31. plot(ecdf(y))
  32. pts=FYs=c(350:700)
  33. for (i in 1:length(pts)) FYs[i]=FY(pts[i],mu*n,sig*sqrt(n))
  34. lines(pts[!is.na(FYs)],FYs[!is.na(FYs)],col="red")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement