# Untitled

Apr 24th, 2021
1,012
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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) {
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")
RAW Paste Data