Wednesday, 20 March 2019

ComplexChooserOption

ComplexChooserOption
bench <- function(N = 1e5, K = 10) {
+     x <- rnorm(N)
+     gc()
+     t <- c()
+     t[1] <- system.time(for (k in 1:K) median(x))[3]
+     t[2] <- system.time(for (k in 1:K) weightedMedian(x))[3]
+     t <- t / t[1]
+     names(t) <- c("median", "weightedMedian")
+     t
+ }
> library(fOptions)






package ‘fExoticOptions’ was built under R version 3.5.1
> library(fExoticOptions)
> FEInDomesticFXOption(TypeFlag = "c", S = 100, E = 1.5,
+                      X = 160, Time = 0.5, r = 0.08, q = 0.05, sigmaS = 0.20,
+                      sigmaE = 0.12, rho = 0.45)

Title:
 FE In Domestic FX Option

Call:
 FEInDomesticFXOption(TypeFlag = "c", S = 100, E = 1.5, X = 160,
     Time = 0.5, r = 0.08, q = 0.05, sigmaS = 0.2, sigmaE = 0.12,
     rho = 0.45)

Parameters:
          Value:
 TypeFlag c   
 S        100 
 E        1.5 
 X        160 
 Time     0.5 
 r        0.08 
 q        0.05 
 sigmaS   0.2 
 sigmaE   0.12 
 rho      0.45 

Option Price:
 8.305561

Description:
 Wed Mar 20 16:24:40 2019

> QuantoOption(TypeFlag = "c", S = 100, Ep = 1.5, X = 105,
+              Time = 0.5, r = 0.08, rf = 0.05, q = 0.04, sigmaS= 0.2,
+              sigmaE = 0.10, rho = 0.30)

Title:
 Quanto Option

Call:
 QuantoOption(TypeFlag = "c", S = 100, Ep = 1.5, X = 105, Time = 0.5,
     r = 0.08, rf = 0.05, q = 0.04, sigmaS = 0.2, sigmaE = 0.1,
     rho = 0.3)

Parameters:
          Value:
 TypeFlag c   
 S        100 
 Ep       1.5 
 X        105 
 Time     0.5 
 r        0.08 
 rf       0.05 
 q        0.04 
 sigmaS   0.2 
 sigmaE   0.1 
 rho      0.3 

Option Price:
 5.328035

Description:
 Wed Mar 20 16:25:02 2019

> EquityLinkedFXOption(TypeFlag = "p", E = 1.5, S = 100,
+                      X = 1.52, Time = 0.25, r = 0.08, rf = 0.05, q = 0.04,
+                      sigmaS = 0.20, sigmaE = 0.12, rho = -0.40)

Title:
 Equity Linked FX Option

Call:
 EquityLinkedFXOption(TypeFlag = "p", E = 1.5, S = 100, X = 1.52,
     Time = 0.25, r = 0.08, rf = 0.05, q = 0.04, sigmaS = 0.2,
     sigmaE = 0.12, rho = -0.4)

Parameters:
          Value:
 TypeFlag p   
 E        1.5 
 S        100 
 X        1.52 
 Time     0.25 
 r        0.08 
 rf       0.05 
 q        0.04 
 sigmaS   0.2 
 sigmaE   0.12 
 rho      -0.4 

Option Price:
 4.208856

Description:
 Wed Mar 20 16:25:20 2019

> TakeoverFXOption(V = 100, B = 100, E = 1.5, X = 1.55, Time = 1,
+                  r = 0.08, rf = 0.06, sigmaV = 0.20, sigmaE = 0.25, rho = 0.1)

Title:
 Takeover FX Option

Call:
 TakeoverFXOption(V = 100, B = 100, E = 1.5, X = 1.55, Time = 1,
     r = 0.08, rf = 0.06, sigmaV = 0.2, sigmaE = 0.25, rho = 0.1)

Parameters:
        Value:                                                                                         
 V      100                                                                                           
 B      100                                                                                           
 E      1.5                                                                                           
 X      1.55                                                                                           
 Time   1                                                                                             
 r      0.08                                                                                           
 rf     0.06                                                                                           
 q      function (save = "default", status = 0, runLast = TRUE) \n.Internal(quit(save, status, runLast))
 sigmaV 0.2                                                                                           
 sigmaE 0.25                                                                                           
 rho    0.1                                                                                           

Option Price:
 4.977936

Description:
 Wed Mar 20 16:25:35 2019

> vanilla <- GBSOption(TypeFlag = "c", S = 100, X = 90, Time = 1,
+                      r = 0.02, b = -0.02, sigma = 0.3)
> KO <- sapply(100:300, FUN = StandardBarrierOption, TypeFlag = "cuo",
+              S = 100, X = 90, K = 0, Time = 1, r = 0.02, b = -0.02, sigma = 0.30)
Warning messages:
1: In if (X >= H) { :
  the condition has length > 1 and only the first element will be used
2: In if (X < H) { :
  the condition has length > 1 and only the first element will be used
> plot(KO[[1]]@price, type = "l",
+      xlab = "barrier distance from spot",
+      ylab = "price of option",
+      main = "Price of KO converges to plain vanilla")
> abline(h = vanilla@price, col = "red")
> library(plot3D)
Warning message:
package ‘plot3D’ was built under R version 3.5.2
> BS_surface <- function(S, Time, FUN, ...) {
+     require(plot3D)
+     n <- length(S)
+     k <- length(Time)
+     m <- matrix(0, n, k)
+     for (i in 1:n){
+         for (j in 1:k){
+             l <- list(S = S[i], Time = Time[j], ...)
+             m[i,j] <- max(do.call(FUN, l)@price, 0)
+         }
+     }
+ persp3D(z = m, xlab = "underlying", ylab = "Remaining time",
+         zlab = "option price", phi = 30, theta = 20, bty = "b2")
+ }
>
> BS_surface(seq(1, 200,length = 200), seq(0, 2, length = 200),
+            GBSOption, TypeFlag = "c", X = 90, r = 0.02, b = 0, sigma = 0.3)
> ExecutiveStockOption(TypeFlag = "c", S = 60, X = 64, Time = 2,
+                      r = 0.07, b = 0.07-0.03, sigma = 0.38, lambda = 0.15)

Title:
 Executive Stock Option Valuation

Call:
 ExecutiveStockOption(TypeFlag = "c", S = 60, X = 64, Time = 2,
     r = 0.07, b = 0.07 - 0.03, sigma = 0.38, lambda = 0.15)

Parameters:
          Value:
 TypeFlag c   
 S        60   
 X        64   
 Time     2   
 r        0.07 
 b        0.04 
 sigma    0.38 
 lambda   0.15 

Option Price:
 9.124388

Description:
 Wed Mar 20 16:51:26 2019

> ForwardStartOption(TypeFlag = "c", S = 60, alpha = 1.1,
+                    time1 = 1, Time2 = 1/4, r = 0.08, b = 0.08-0.04, sigma = 0.30)

Title:
 Forward Start Option Valuation

Call:
 ForwardStartOption(TypeFlag = "c", S = 60, alpha = 1.1, time1 = 1,
     Time2 = 1/4, r = 0.08, b = 0.08 - 0.04, sigma = 0.3)

Parameters:
          Value:
 TypeFlag c   
 S        60   
 alpha    1.1 
 time1    1   
 Time2    0.25 
 r        0.08 
 b        0.04 
 sigma    0.3 

Option Price:
 4.406449

Description:
 Wed Mar 20 16:52:06 2019

> RatchetOption(TypeFlag = "c", S = 60, alpha = 1.1, time1 = c(1.00, 0.75),
+               Time2 = c(0.75, 0.50), r = 0.08, b = 0.04, sigma = 0.30)

Title:
 Ratchet Option Valuation

Call:
 RatchetOption(TypeFlag = "c", S = 60, alpha = 1.1, time1 = c(1,
     0.75), Time2 = c(0.75, 0.5), r = 0.08, b = 0.04, sigma = 0.3)

Parameters:
          Value:
 TypeFlag c   
 S        60   
 alpha    1.1 
 time11   1   
 time12   0.75 
 Time21   0.75 
 Time22   0.5 
 r        0.08 
 b        0.04 
 sigma    0.3 

Option Price:
 3.2132

Description:
 Wed Mar 20 16:52:51 2019

> TimeSwitchOption(TypeFlag = "c", S = 100, X = 110, Time = 1,
+                  r = 0.06, b = 0.06, sigma = 0.26, A = 5, m = 0, dt = 1/365)

Title:
 Time Switch Option Valuation

Call:
 TimeSwitchOption(TypeFlag = "c", S = 100, X = 110, Time = 1,
     r = 0.06, b = 0.06, sigma = 0.26, A = 5, m = 0, dt = 1/365)

Parameters:
          Value:           
 TypeFlag c                 
 S        100               
 X        110               
 Time     1                 
 r        0.06             
 b        0.06             
 sigma    0.26             
 A        5                 
 m        0                 
 d        0.00273972602739726

Option Price:
 1.375037

Description:
 Wed Mar 20 16:53:19 2019

> SimpleChooserOption(S = 50, X = 50, time1 = 1/4, Time2 = 1/2,
+                     r = 0.08, b = 0.08, sigma = 0.25)

Title:
 Simple Chooser Option Valuation

Call:
 SimpleChooserOption(S = 50, X = 50, time1 = 1/4, Time2 = 1/2,
     r = 0.08, b = 0.08, sigma = 0.25)

Parameters:
       Value:
 S     50   
 X     50   
 time1 0.25 
 Time2 0.5 
 r     0.08 
 b     0.08 
 sigma 0.25 

Option Price:
 6.107073

Description:
 Wed Mar 20 16:53:35 2019

> ComplexChooserOption(S = 50, Xc = 55, Xp = 48, Time = 0.25,
+                      Timec = 0.50, Timep = 0.5833, r = 0.10, b = 0.1-0.05,
+                      sigma = 0.35, doprint = TRUE)

Critical Value:
[1] 51.11561


Title:
 Complex Chooser Option Valuation

Call:
 ComplexChooserOption(S = 50, Xc = 55, Xp = 48, Time = 0.25, Timec = 0.5,
     Timep = 0.5833, r = 0.1, b = 0.1 - 0.05, sigma = 0.35, doprint = TRUE)

Parameters:
               Value:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
 S             50                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 Xc            55                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 Xp            48                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 time          new("standardGeneric", .Data = function (x, ...) \nstandardGeneric("time"), generic = "time", package = "stats", group = list(), valueClass = character(0), signature = "x", default = new("derivedDefaultMethod", .Data = function (x, ...) \nUseMethod("time"), target = new("signature", .Data = "ANY", names = "x", package = "methods"), defined = new("signature", .Data = "ANY", names = "x", package = "methods"), generic = "time"), skeleton = (new("derivedDefaultMethod", .Data = function (x, ...) \nUseMethod("time"), target = new("signature", .Data = "ANY", names = "x", package = "methods"), defined = new("signature", .Data = "ANY", names = "x", package = "methods"), generic = "time"))(x, ...))
 Timec         0.5                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 Timep         0.5833                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
 r             0.1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 b             0.05                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
 sigma         0.35                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
 criticalValue 51.1156056231326                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       

Option Price:
 6.050727

Description:
 Wed Mar 20 16:53:59 2019

> OptionOnOption(TypeFlag = "pc", S = 500, X1 = 520, X2 = 50,
+                time1 = 1/2, Time2 = 1/4, r = 0.08, b = 0.08-0.03, sigma = 0.35)

Title:
 Option On Option Valuation

Call:
 OptionOnOption(TypeFlag = "pc", S = 500, X1 = 520, X2 = 50, time1 = 1/2,
     Time2 = 1/4, r = 0.08, b = 0.08 - 0.03, sigma = 0.35)

Parameters:
               Value:         
 S             500           
 X1            520           
 X2            50             
 time1         0.5           
 Time2         0.25           
 r             0.08           
 b             0.05           
 sigma         0.35           
 criticalValue 538.316546560699

Option Price:
 21.19647

Description:
 Wed Mar 20 16:54:18 2019

> HolderExtendibleOption(TypeFlag = "c", S = 100, X1 = 100,
+                        X2 = 105, time1 = 0.50, Time2 = 0.75, r = 0.08, b = 0.08,
+                        sigma = 0.25, A = 1)

Title:
 Holder Extendible Option Valuation

Call:
 HolderExtendibleOption(TypeFlag = "c", S = 100, X1 = 100, X2 = 105,
     time1 = 0.5, Time2 = 0.75, r = 0.08, b = 0.08, sigma = 0.25,
     A = 1)

Parameters:
          Value:
 TypeFlag c   
 S        100 
 X1       100 
 X2       105 
 time1    0.5 
 Time2    0.75 
 r        0.08 
 b        0.08 
 sigma    0.25 
 A        1   

Option Price:
 2.729216

Description:
 Wed Mar 20 16:54:38 2019

> WriterExtendibleOption(TypeFlag = "c", S = 80, X1 = 90, X2 = 82,
+                        time1 = 0.50, Time2 = 0.75, r = 0.10, b = 0.10, sigma = 0.30)

Title:
 Writer Extendible Option Valuation

Call:
 WriterExtendibleOption(TypeFlag = "c", S = 80, X1 = 90, X2 = 82,
     time1 = 0.5, Time2 = 0.75, r = 0.1, b = 0.1, sigma = 0.3)

Parameters:
          Value:
 TypeFlag c   
 S        80   
 X1       90   
 X2       82   
 time1    0.5 
 Time2    0.75 
 r        0.1 
 b        0.1 
 sigma    0.3 

Option Price:
 6.823755

Description:
 Wed Mar 20 16:54:55 2019

> install.packages("plot3Drgl")
Installing package into ‘C:/Users/ADMIN/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/plot3Drgl_1.0.1.zip'
Content type 'application/zip' length 248378 bytes (242 KB)
downloaded 242 KB

package ‘plot3Drgl’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\ADMIN\AppData\Local\Temp\RtmpySBgIa\downloaded_packages
> library(plot3D-packages)
Error in library(plot3D - packages) : 'package' must be of length 1
> par(mfrow = c(2, 1))
> x <- rchisq(1000, df = 4)
> hs <- hist(x, breaks = 15)
> hist3D(x = hs$mids, y = 1, z = matrix(ncol = 1, data = hs$density),
+        bty = "g", ylim = c(0., 2.0), scale = FALSE, expand = 20,
+        border = "black", col = "white", shade = 0.3, space = 0.1,
+        theta = 20, phi = 20, main = "3-D perspective")
> plotdev(xlim = c(-0.2, 1.2), ylim = c(-0.2, 1.2), theta = 45)
> ribbon3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
+          col = "lightblue", border = "black", shade = 0.2, ltheta = 90,
+          space = 0.3, ticktype = "detailed", d = 2, curtain = TRUE)
Error in nrow(z) : object 'VV' not found
> pm <- par("mfrow")
>
> par(mfrow = c(2, 2))
> persp3D(z = volcano, main = "volcano", clab = c("height", "m"),
+         breaks = seq(80,200, by = 10))
> window()
Error in hasTsp(x) : argument "x" is missing, with no default
> persp3D(z = volcano, main = "volcano", clab = c("height", "m"),
+         breaks = seq(80,200, by = 10))
> persp3D(z = volcano, x = 1: nrow(volcano), y = 1:ncol(volcano),
+         expand = 0.3, main = "volcano", facets = FALSE, scale = FALSE,
+         clab = "height, m", colkey = list(side = 1, length = 0.5))
> V <- volcano[, seq(1, ncol(volcano), by = 3)]  # lower resolution
> ribbon3D(z = V, colkey = list(width = 0.5, length = 0.5,
+                               cex.axis = 0.8, side = 2), clab = "m")
> Vy <- volcano[seq(1, nrow(volcano), by = 3), ]
> ribbon3D(z = Vy, expand = 0.3, space = 0.3, along = "y",
+          colkey = list(width = 0.5, length = 0.5, cex.axis = 0.8))
>
> x <- seq(-pi, pi, by = 0.2)
> y <- seq(-pi, pi, by = 0.3)
> grid <- mesh(x, y)
> z    <- with(grid, cos(x) * sin(y))
> par(mfrow = c(2,2))
> persp3D(z = z, x = x, y = y)
> persp3D(z = z, x = x, y = y, facets = FALSE, curtain = TRUE)
> ribbon3D(z = z, x = x, y = y, along = "xy", space = 0.3)
> hist3D(z = z, x = x, y = y, border = "black")
> par(mfrow = c(2, 2))
> x <- seq(1, nrow(volcano), by = 3)
> y <- seq(1, ncol(volcano), by = 3)
> Volcano <- volcano [x, y]
> ribbon3D(z = Volcano, contour = TRUE, zlim= c(-100, 200),
+          image = TRUE)
> persp3D(z = Volcano, contour = TRUE, zlim= c(-200, 200), image = FALSE)
> persp3D(z = Volcano, x = x, y = y, scale = FALSE,
+         contour = list(nlevels = 20, col = "red"),
+         zlim = c(-200, 200), expand = 0.2,
+         image = list(col = grey (seq(0, 1, length.out = 100))))
> persp3D(z = Volcano, contour = list(side = c("zmin", "z", "350")),
+         zlim = c(-100, 400), phi = 20, image = list(side = 350))
> par(mfrow = c(2, 2))
> persp3D(z = Volcano, shade = 0.5, colkey = FALSE)
> persp3D(z = Volcano, inttype = 2, shade = 0.5, colkey = FALSE)
> x <- y <- seq(0, 2*pi, length.out = 10)
> z <- with (mesh(x, y), cos(x) *sin(y)) + runif(100)
> cv <- matrix(nrow = 10, 0.5*runif(100))
> persp3D(x, y, z, colvar = cv)
> persp3D(x, y, z, colvar = cv, inttype = 2)
> par(mfrow = c(2, 2))
> VV <- V2 <- volcano[10:15, 10:15]
> V2[3:4, 3:4] <- NA
> V2[4, 5] <- NA
> image2D(V2, border = "black")
> persp3D(z = VV, colvar = V2, inttype = 1, theta = 0,
+         phi = 20, border = "black", main = "inttype = 1")
> persp3D(z = VV, colvar = V2, inttype = 2, theta = 0,
+         phi = 20, border = "black", main = "inttype = 2")
> persp3D(z = VV, colvar = V2, inttype = 3, theta = 0,
+         phi = 20, border = "black", main = "inttype = 3")
> par(mfrow = c(1, 1))
> panelfirst <- function(trans) {
+     zticks <- seq(100, 180, by = 20)
+     len <- length(zticks)
+     XY0 <- trans3D(x = rep(1, len), y = rep(1, len), z = zticks,
+                    pmat = trans)
+     XY1 <- trans3D(x = rep(1, len), y = rep(61, len), z = zticks,
+                    pmat = trans)
+     segments(XY0$x, XY0$y, XY1$x, XY1$y, lty = 2)
+   
+     rm <- rowMeans(volcano)           
+     XY <- trans3D(x = 1:87, y = rep(ncol(volcano), 87),
+                   z = rm, pmat = trans)
+     lines(XY, col = "blue", lwd = 2)
+ }
> persp3D(z = volcano, x = 1:87, y = 1: 61, scale = FALSE, theta = 10,
+         expand = 0.2, panel.first = panelfirst, colkey = FALSE)
> par(mfrow = c(2, 2))
> persp3D(z = volcano, shade = 0.3, col = gg.col(100))
> persp3D(z = volcano, lighting = TRUE, lphi = 90)
> persp3D(z = volcano, col = "lightblue", colvar = NULL,
+         shade = 0.3, bty = "b2")
> persp3D(z = volcano, col = "lightblue", colvar = NULL,
+         shade = 0.3, bty = "b2")
> volcx <- matrix(nrow = 87, ncol = 61, data = 1:87)
> volcx <- volcx + matrix(nrow = 87, ncol = 61,
+                         byrow = TRUE, data = seq(0., 15, length.out = 61))
> volcy <- matrix(ncol = 87, nrow = 61, data = 1:61)
> volcy <- t(volcy + matrix(ncol = 87, nrow = 61,
+                           byrow = TRUE, data = seq(0., 15, length.out = 87)))
> persp3D(volcano, x = volcx, y = volcy, phi = 80)
> par(mfrow = c(1, 1))
> clim <- range(volcano)
> persp3D(z = volcano, zlim = c(100, 600), clim = clim,
+         box = FALSE, plot = FALSE)
> persp3D(z = volcano, zlim = c(100, 600), clim = clim,
+         box = FALSE, plot = FALSE)
> persp3D(z = volcano + 200, clim = clim, colvar = volcano,
+         add = TRUE, colkey = FALSE, plot = FALSE)
> persp3D(z = volcano + 400, clim = clim, colvar = volcano,
+         add = TRUE, colkey = FALSE)
> par(mfrow = c(2, 2))
> VV   <- volcano[seq(1, 87, 15), seq(1, 61, 15)]
> hist3D(z = VV, scale = FALSE, expand = 0.01, border = "black")
> hist3D(z = VV, scale = FALSE, expand = 0.01,
+        alpha = 0.5, opaque.top = TRUE, border = "black")
> hist3D(z = VV, scale = FALSE, expand = 0.01, facets = FALSE, lwd = 2)
> hist3D(z = VV, scale = FALSE, expand = 0.01, facets = NA)
> par(mfrow = c(2, 2))
> hist3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
+        col = "#0072B2", border = "black", shade = 0.2, ltheta = 90,
+        space = 0.3, ticktype = "detailed", d = 2)
> plotdev(xlim = c(-0.2, 1.2), ylim = c(-0.2, 1.2), theta = 45)
>
> ribbon3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
+          col = "lightblue", border = "black", shade = 0.2, ltheta = 90,
+          space = 0.3, ticktype = "detailed", d = 2, curtain = TRUE)
> ribbon3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20, zlim = c(95,183),
+          col = "lightblue", lighting = TRUE, ltheta = 50, along = "y",
+          space = 0.7, ticktype = "detailed", d = 2, curtain = TRUE)
>
> par(mfrow = c(2, 1))
> x <- rchisq(1000, df = 4)
> hs <- hist(x, breaks = 15)
> hist3D(x = hs$mids, y = 1, z = matrix(ncol = 1, data = hs$density),
+        bty = "g", ylim = c(0., 2.0), scale = FALSE, expand = 20,
+        border = "black", col = "white", shade = 0.3, space = 0.1,
+        theta = 20, phi = 20, main = "3-D perspective")
> par(mfrow = pm)
> pm <- par("mfrow")
> pmar <- par("mar")
> par(mfrow = c(2, 2))
> x <- y <- z <- seq(-1, 1, by = 0.1)
> grid   <- mesh(x, y, z)
> colvar <- with(grid, x*exp(-x^2 - y^2 - z^2))
> slice3D  (x, y, z, colvar = colvar, theta = 60)
> slicecont3D (x, y, z, ys = seq(-1, 1, by = 0.5), colvar = colvar,
+              theta = 60, border = "black")
> slice3D  (x, y, z, xs = c(-1, -0.5, 0.5), ys = c(-1, 0, 1),
+           zs = c(-1, 0), colvar = colvar,
+           theta = 60, phi = 40)
> XY <- mesh(x, y)
> ZZ <- XY$x*XY$y
> slice3D  (x, y, z, xs = XY$x, ys = XY$y, zs = ZZ, colvar = colvar,
+           lighting =  TRUE, lphi = 90, ltheta = 0)
> par(mfrow = c(1, 1))
> x <- y <- z <- seq(-4, 4, by = 0.2)
> M <- mesh(x, y, z)
> R <- with (M, sqrt(x^2 + y^2 + z^2))
> p <- sin(2*R) /(R+1e-3)
> slice3D(x, y, z, colvar = p, col = jet.col(alpha = 0.5),
+         xs = 0, ys = c(-4, 0, 4), zs = NULL, d = 2)
> slice3D(x, y, z, colvar = p, d = 2, theta = 60, border = "black",
+         xs = c(-4, 0), ys = c(-4, 0, 4), zs = c(-4, 0))
> slice3D(x, y, z, colvar = p, d = 2, theta = 60, border = "green",
+         xs = c(-4, 0), ys = c(-4, 0, 4), zs = c(-4, 0))
> slice3D(x, y, z, colvar = p, d = 2, theta = 60, border = "red",
+         xs = c(-4, 0), ys = c(-4, 0, 4), zs = c(-4, 0))
> data(Oxsat)
> Ox <- Oxsat$val[,  Oxsat$lat > - 5 & Oxsat$lat < 5, ]
> slice3D(x = Oxsat$lon, z = -Oxsat$depth, y = 1:5, colvar = Ox,
+         ys = 1:5, zs = NULL, NAcol = "black",
+         expand = 0.4, theta = 45, phi = 45)
> par(mfrow = c(2, 2), mar  = c(2, 2, 2, 2))
> x <- y <- z <- seq(-2, 2, length.out = 15)
> F <- with(xyz, log(x^2 + y^2 + z^2 +
+                        10*(x^2 + y^2) * (y^2 + z^2) ^2))
Error in with(xyz, log(x^2 + y^2 + z^2 + 10 * (x^2 + y^2) * (y^2 + z^2)^2)) :
  object 'xyz' not found
> xyz <- mesh(x, y, z)
> F <- with(xyz, log(x^2 + y^2 + z^2 +
+                        10*(x^2 + y^2) * (y^2 + z^2) ^2))
> isosurf3D(x, y, z, F, level = 1, shade = 0.9,
+           col = "yellow", border = "orange")
> isosurf3D(x, y, z, F, level = 2, lighting = TRUE,
+           lphi = 0, ltheta = 0, col = "blue", shade = NA)
> isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
+           col = c("red", "blue", "yellow"),
+           clab = "F", alpha = 0.2, theta = 0, lighting = TRUE)
> isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
+           col = c("red", "blue", "yellow"),
+           shade = NA, plot = FALSE, clab = "F") 
> plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
> iso <- createisosurf(x, y, z, F, level = 2)
> head(iso)
              x          y         z
[1,]  0.0000000 -0.1179802 -2.000000
[2,] -0.1204879  0.0000000 -2.000000
[3,]  0.0000000 -0.2073773 -1.714286
[4,]  0.0000000 -0.2073773 -1.714286
[5,] -0.1204879  0.0000000 -2.000000
[6,] -0.2138894  0.0000000 -1.714286
> triangle3D(iso, col = "green", shade = 0.3)
> x <- y <- z <- seq(-2, 2, length.out = 50)
> xyz <- mesh(x, y, z)
> F <- with(xyz, log(x^2 + y^2 + z^2 +
+                        10*(x^2 + y^2) * (y^2 + z^2) ^2))
> isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
+           col = c("red", "blue", "yellow"),
+           shade = NA, plot = FALSE, clab = "F") 
> plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
> par(mfrow = c(2, 2), mar  = c(2, 2, 2, 2))
> x <- y <- z <- seq(-2, 2, length.out = 70)
> xyz <- mesh(x, y, z)
> F <- with(xyz, log(x^2 + y^2 + z^2 +
+                        10*(x^2 + y^2) * (y^2 + z^2) ^2))
> voxel3D(x, y, z, F, level = 4, pch = ".", cex = 5)
> plotdev(theta = 45, phi = 0)
> plotdev(theta = 30, phi = 0)
> plotdev(theta = 90, phi = 10)
> plotdev(theta = 90, phi = 25)
> plotdev(theta = 120, phi = 20)
> vox <- createvoxel(x, y, z, F, level = 4)
> scatter3D(vox$x, vox$y, vox$z, colvar = vox$y,
+           bty = "g", colkey = FALSE)
>
> par(mfrow = c(1, 1), mar = c(2, 2, 2, 2))
>
> Hypox <- createvoxel(Oxsat$lon, Oxsat$lat, Oxsat$depth[1:19],
+                      Oxsat$val[,,1:19], level = 40, operator = "<")
> panel <- function(pmat) {  # an image at the bottom
+     Nx <- length(Oxsat$lon)
+     Ny <- length(Oxsat$lat)
+     M <- mesh(Oxsat$lon, Oxsat$lat)
+     xy <- trans3D(pmat = pmat, x = as.vector(M$x), y = as.vector(M$y),
+                   z = rep(-1000, length.out = Nx*Ny))
+     x <- matrix(nrow = Nx, ncol = Ny, data = xy$x)
+     y <- matrix(nrow = Nx, ncol = Ny, data = xy$y)
+     Bat <- Oxsat$val[,,1]; Bat[!is.na(Bat)] <- 1
+     image2D(x = x, y = y, z = Bat, NAcol = "black", col = "grey",
+             add = TRUE, colkey = FALSE)
+ }
> scatter3D(Hypox$x, Hypox$y, -Hypox$z, colvar = Hypox$cv,
+           panel.first = panel, pch = ".", bty = "b",
+           theta = 30, phi = 20, ticktype = "detailed",
+           zlim = c(-1000,0), xlim = range(Oxsat$lon),
+           ylim = range(Oxsat$lat) )
> pm <- par("mfrow")
> par(mfrow = c(1, 1))
> image2D(z = Oxsat$val[ , , 1], x = Oxsat$lon, y = Oxsat$lat,
+         main = "surface oxygen saturation (%) for 2005")
> image2D(z = Oxsat$val[ , , 1], x = Oxsat$lon, y = Oxsat$lat,
+         main = "surface oxygen saturation (%) for 2018")
> lon <- Oxsat$lon
> image2D (z = Oxsat$val, margin = c(2, 3), x = Oxsat$lat,
+          y = Oxsat$depth, subset = (lon > 18 & lon < 23),
+          ylim = c(5500, 0), NAcol = "black", zlim = c(0, 110),
+          xlab = "latitude", ylab = "depth, m")
> ImageOcean()
> abline ( v = lon[lon > 18 & lon < 23])
> par(mfrow = c(1, 1))
> ii <- which (Oxsat$lon > -90 & Oxsat$lon < 90)
> jj <- which (Oxsat$lat > 0 & Oxsat$lat < 90)
> xs <- Oxsat$lon[ii[length(ii)]]
> ys <- Oxsat$lat[jj[1]]
> slice3D(colvar = Oxsat$val[ii,jj,], x = Oxsat$lon[ii], 
+         y = Oxsat$lat[jj], z = -Oxsat$depth,
+         NAcol = "black", xs = xs, ys = ys, zs = 0,
+         theta = 35, phi = 50, colkey = list(length = 0.5),
+         expand = 0.5, ticktype = "detailed",
+         clab = "%", main = "Oxygen saturation",
+         xlab = "longitude", ylab = "latitude", zlab = "depth")
> x <- y <- z <- c(-1 , 0, 1)
> (M <- mesh(x, y, z))
$`x`
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 2

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 3

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1


$y
, , 1

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 2

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 3

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1


$z
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]   -1   -1   -1
[3,]   -1   -1   -1

, , 2

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

, , 3

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1


> V <- with (M, x/2 * sin(y) *sqrt(z+2))
> scatter3D(M$x, M$y, M$z, V, pch = "+", cex = 3, colkey = FALSE)
> x <- runif(10)
> y <- runif(10)
> polygon2D(x, y)
> x <- runif(59)
> y <- runif(59)
> ii <- seq(5, 55, by  = 5)
> x[ii] <- y[ii] <- NA
> polygon2D(x, y, border = "black", lwd = 3,
+           colvar = 1:(length(ii) + 1),
+           col = gg.col(), alpha = 0.2,
+           main = "polygon2D")
> par(mfrow = c(2, 2))
> x <- y <- z <- seq(-1, 1, by = 0.1)
> grid   <- mesh(x, y, z)
> colvar <- with(grid, x*exp(-x^2 - y^2 - z^2))
> slice3D  (x, y, z, colvar = colvar, theta = 60)
> slicecont3D (x, y, z, ys = seq(-1, 1, by = 0.5), colvar = colvar,
+              theta = 60, border = "black")
> slice3D  (x, y, z, xs = c(-1, -0.5, 0.5), ys = c(-1, 0, 1),
+           zs = c(-1, 0), colvar = colvar,
+           theta = 60, phi = 40)
> slice3D  (x, y, z, xs = c(-1, -0.5, 0.5), ys = c(-1, 0, 1),
+           zs = c(-1, 0), colvar = colvar,
+           theta = 90, phi = 70)

Monday, 18 March 2019

StandardBarrierOption

StandardBarrierOption
[ reached getOption("max.print") -- omitted 831 rows ]

$R
            [,1]     [,2]     [,3]     [,4]
   [1,] 3.524952 3.476878 3.409521 3.290123
   [2,] 3.514023 3.466389 3.399649 3.281343
   [3,] 3.527616 3.479434 3.411928 3.292263
   [4,] 3.516217 3.468495 3.401631 3.283105
   [5,] 3.527133 3.478971 3.411492 3.291876
   [6,] 3.512964 3.465372 3.398692 3.280492
   [7,] 3.480904 3.434603 3.369731 3.254733
   [8,] 3.445434 3.400560 3.337688 3.226234
   [9,] 3.452784 3.407615 3.344328 3.232140
  [10,] 3.470931 3.425032 3.360721 3.246721
  [11,] 3.439446 3.394814 3.332279 3.221424
  [12,] 3.449771 3.404723 3.341606 3.229719
  [13,] 3.483697 3.437283 3.372253 3.256977
  [14,] 3.494716 3.447859 3.382208 3.265830
  [15,] 3.521277 3.473350 3.406201 3.287170
  [16,] 3.495879 3.448976 3.383258 3.266765
  [17,] 3.491624 3.444891 3.379414 3.263346
  [18,] 3.509198 3.461758 3.395290 3.277466
  [19,] 3.515956 3.468244 3.401395 3.282896
  [20,] 3.509752 3.462290 3.395790 3.277911
  [21,] 3.528847 3.480616 3.413040 3.293252
  [22,] 3.534825 3.486354 3.418440 3.298056
  [23,] 3.519893 3.472023 3.404951 3.286059
  [24,] 3.481423 3.435101 3.370199 3.255150
  [25,] 3.491160 3.444446 3.378995 3.262973
  [26,] 3.491862 3.445120 3.379629 3.263537
  [27,] 3.492130 3.445377 3.379871 3.263753
  [28,] 3.480309 3.434032 3.369193 3.254255
  [29,] 3.508341 3.460936 3.394516 3.276777
  [30,] 3.503402 3.456196 3.390054 3.272809
  [31,] 3.481586 3.435258 3.370346 3.255281
  [32,] 3.478578 3.432371 3.367629 3.252864
  [33,] 3.466595 3.420870 3.356804 3.243236
  [34,] 3.483430 3.437028 3.372012 3.256763
  [35,] 3.470891 3.424993 3.360685 3.246688
  [36,] 3.452084 3.406943 3.343695 3.231578
  [37,] 3.461109 3.415604 3.351848 3.238828
  [38,] 3.490576 3.443886 3.378467 3.262504
  [39,] 3.510494 3.463001 3.396460 3.278507
  [40,] 3.535382 3.486888 3.418943 3.298503
  [41,] 3.543177 3.494369 3.425985 3.304766
  [42,] 3.551851 3.502694 3.433821 3.311735
  [43,] 3.551489 3.502347 3.433494 3.311444
  [44,] 3.541422 3.492685 3.424400 3.303356
  [45,] 3.579011 3.528761 3.458356 3.333557
  [46,] 3.569951 3.520066 3.450172 3.326278
  [47,] 3.587000 3.536428 3.465573 3.339976
  [48,] 3.611759 3.560190 3.487939 3.359868
  [49,] 3.597672 3.546671 3.475214 3.348550
  [50,] 3.611905 3.560330 3.488071 3.359985
  [51,] 3.584484 3.534014 3.463300 3.337954
  [52,] 3.569083 3.519232 3.449387 3.325580
  [53,] 3.563939 3.514295 3.444740 3.321447
  [54,] 3.572158 3.522184 3.452165 3.328051
  [55,] 3.570999 3.521071 3.451118 3.327120
  [56,] 3.545592 3.496687 3.428167 3.306707
  [57,] 3.570521 3.520612 3.450686 3.326735
  [58,] 3.538524 3.489904 3.421782 3.301028
  [59,] 3.537333 3.488760 3.420706 3.300071
  [60,] 3.544279 3.495427 3.426980 3.305651
  [61,] 3.555120 3.505831 3.436774 3.314362
  [62,] 3.537418 3.488842 3.420782 3.300139
  [63,] 3.556348 3.507010 3.437883 3.315348
  [64,] 3.562840 3.513240 3.443747 3.320564
  [65,] 3.568243 3.518426 3.448629 3.324906
  [66,] 3.578486 3.528256 3.457881 3.333135
  [67,] 3.553296 3.504080 3.435126 3.312896
  [68,] 3.553544 3.504318 3.435350 3.313095
  [69,] 3.551680 3.502530 3.433667 3.311598
  [70,] 3.565942 3.516218 3.446550 3.323057
  [71,] 3.552652 3.503462 3.434544 3.312379
  [72,] 3.536197 3.487670 3.419680 3.299158
  [73,] 3.554588 3.505321 3.436294 3.313935
  [74,] 3.538448 3.489830 3.421713 3.300967
  [75,] 3.566568 3.516819 3.447116 3.323560
  [76,] 3.566252 3.516516 3.446831 3.323306
  [77,] 3.551813 3.502658 3.433787 3.311705
  [78,] 3.581978 3.531608 3.461036 3.335941
  [79,] 3.567565 3.517775 3.448016 3.324361
  [80,] 3.585303 3.534800 3.464040 3.338612
  [81,] 3.602158 3.550976 3.479266 3.352154
  [82,] 3.595775 3.544849 3.473500 3.347025
  [83,] 3.621166 3.569219 3.496437 3.367426
  [84,] 3.612526 3.560927 3.488633 3.360485
  [85,] 3.630367 3.578049 3.504749 3.374819
  [86,] 3.640723 3.587988 3.514104 3.383139
  [87,] 3.629787 3.577492 3.504225 3.374352
  [88,] 3.597507 3.546513 3.475065 3.348418
  [89,] 3.616671 3.564904 3.492376 3.363814
  [90,] 3.645939 3.592994 3.518816 3.387330
  [91,] 3.640392 3.587671 3.513805 3.382873
  [92,] 3.646505 3.593538 3.519327 3.387785
  [93,] 3.649038 3.595969 3.521616 3.389820
  [94,] 3.642831 3.590012 3.516008 3.384833
  [95,] 3.613739 3.562091 3.489728 3.361459
  [96,] 3.616917 3.565141 3.492599 3.364012
  [97,] 3.626998 3.574815 3.501705 3.372111
  [98,] 3.641588 3.588819 3.514886 3.383834
  [99,] 3.606377 3.555025 3.483077 3.355544
 [100,] 3.584083 3.533629 3.462938 3.337632
 [101,] 3.564228 3.514573 3.445002 3.321679
 [102,] 3.569536 3.519667 3.449796 3.325944
 [103,] 3.549071 3.500026 3.431309 3.309501
 [104,] 3.514779 3.467114 3.400332 3.281950
 [105,] 3.509099 3.461663 3.395201 3.277387
 [106,] 3.496934 3.449988 3.384211 3.267612
 [107,] 3.502696 3.455518 3.389417 3.272242
 [108,] 3.521565 3.473627 3.406462 3.287402
 [109,] 3.545473 3.496572 3.428059 3.306611
 [110,] 3.534466 3.486009 3.418116 3.297767
 [111,] 3.526632 3.478491 3.411039 3.291473
 [112,] 3.524769 3.476702 3.409356 3.289976
 [113,] 3.534014 3.485575 3.417707 3.297404
 [114,] 3.561572 3.512024 3.442603 3.319546
 [115,] 3.553364 3.504147 3.435188 3.312951
 [116,] 3.557229 3.507856 3.438679 3.316056
 [117,] 3.535744 3.487236 3.419271 3.298794
 [118,] 3.498832 3.451809 3.385926 3.269137
 [119,] 3.505033 3.457761 3.391528 3.274120
 [120,] 3.493463 3.446657 3.381076 3.264824
 [121,] 3.466470 3.420750 3.356691 3.243136
 [122,] 3.470547 3.424663 3.360374 3.246412
 [123,] 3.484611 3.438161 3.373079 3.257711
 [124,] 3.460320 3.414848 3.351136 3.238195
 [125,] 3.458877 3.413462 3.349832 3.237035
 [126,] 3.470041 3.424178 3.359917 3.246005
 [127,] 3.474500 3.428457 3.363945 3.249588
 [128,] 3.498821 3.451798 3.385915 3.269128
 [129,] 3.499189 3.452152 3.386248 3.269424
 [130,] 3.479347 3.433109 3.368324 3.253482
 [131,] 3.463421 3.417823 3.353937 3.240686
 [132,] 3.443415 3.398623 3.335864 3.224613
 [133,] 3.432512 3.388159 3.326015 3.215852
 [134,] 3.451741 3.406613 3.343385 3.231302
 [135,] 3.476802 3.430666 3.366025 3.251437
 [136,] 3.473425 3.427425 3.362974 3.248724
 [137,] 3.464477 3.418837 3.354891 3.241535
 [138,] 3.490583 3.443893 3.378474 3.262510
 [139,] 3.488932 3.442307 3.376982 3.261183
 [140,] 3.477156 3.431006 3.366345 3.251722
 [141,] 3.469077 3.423252 3.359046 3.245231
 [142,] 3.437310 3.392764 3.330350 3.219708
 [143,] 3.430656 3.386378 3.324338 3.214361
 [144,] 3.502938 3.455750 3.389635 3.272436
 [145,] 3.510158 3.462679 3.396157 3.278237
 [146,] 3.500796 3.453694 3.387700 3.270715
 [147,] 3.504509 3.457258 3.391054 3.273699
 [148,] 3.461383 3.415868 3.352096 3.239049
 [149,] 3.465483 3.419803 3.355800 3.242343
 [150,] 3.455177 3.409911 3.346489 3.234063
 [151,] 3.478722 3.432509 3.367759 3.252980
 [152,] 3.454850 3.409597 3.346194 3.233800
 [153,] 3.451062 3.405962 3.342772 3.230756
 [154,] 3.450022 3.404964 3.341833 3.229921
 [155,] 3.440712 3.396029 3.333423 3.222441
 [156,] 3.453443 3.408248 3.344923 3.232670
 [157,] 3.465317 3.419644 3.355650 3.242210
 [158,] 3.449591 3.404551 3.341444 3.229575
 [159,] 3.454714 3.409467 3.346071 3.233691
 [160,] 3.458646 3.413241 3.349624 3.236850
 [161,] 3.458193 3.412806 3.349214 3.236486
 [162,] 3.433115 3.388738 3.326560 3.216337
 [163,] 3.443312 3.398524 3.335771 3.224530
 [164,] 3.407489 3.364144 3.303410 3.195748
 [165,] 3.414723 3.371086 3.309945 3.201560
 [166,] 3.401421 3.358319 3.297928 3.190872
 [167,] 3.376756 3.334647 3.275647 3.171056
 [168,] 3.386623 3.344117 3.284561 3.178983
 [169,] 3.339431 3.298825 3.241930 3.141067
 [170,] 3.361588 3.320090 3.261945 3.158869
 [171,] 3.359966 3.318534 3.260480 3.157566
 [172,] 3.366152 3.324470 3.266068 3.162536
 [173,] 3.333094 3.292743 3.236205 3.135976
 [174,] 3.348537 3.307565 3.250156 3.148383
 [175,] 3.320825 3.280968 3.225122 3.126118
 [176,] 3.305508 3.266267 3.211284 3.113811
 [177,] 3.275905 3.237857 3.184543 3.090027
 [178,] 3.269569 3.231775 3.178819 3.084936
 [179,] 3.293465 3.254709 3.200406 3.104136
 [180,] 3.288308 3.249760 3.195747 3.099992
 [181,] 3.289148 3.250566 3.196506 3.100667
 [182,] 3.325168 3.285136 3.229045 3.129608
 [183,] 3.283788 3.245421 3.191663 3.096360
 [184,] 3.264928 3.227321 3.174627 3.081208
 [185,] 3.260036 3.222626 3.170207 3.077277
 [186,] 3.265602 3.227968 3.175236 3.081749
 [187,] 3.275191 3.237171 3.183898 3.089454
 [188,] 3.263742 3.226183 3.173555 3.080255
 [189,] 3.262178 3.224682 3.172143 3.078999
 [190,] 3.260451 3.223024 3.170582 3.077611
 [191,] 3.236360 3.199903 3.148820 3.058255
 [192,] 3.246054 3.209207 3.157577 3.066044
 [193,] 3.259388 3.222005 3.169622 3.076757
 [194,] 3.225193 3.189186 3.138732 3.049283
 [195,] 3.206821 3.171553 3.122135 3.034522
 [196,] 3.222493 3.186594 3.136293 3.047113
 [197,] 3.205387 3.170177 3.120839 3.033369
 [198,] 3.194897 3.160109 3.111363 3.024941
 [199,] 3.238198 3.201668 3.150480 3.059732
 [200,] 3.236590 3.200124 3.149027 3.058440
 [201,] 3.245985 3.209141 3.157514 3.065988
 [202,] 3.227145 3.191060 3.140495 3.050851
 [203,] 3.229770 3.193578 3.142866 3.052960
 [204,] 3.226351 3.190297 3.139778 3.050213
 [205,] 3.251475 3.214410 3.162474 3.070399
 [206,] 3.284988 3.246574 3.192748 3.097325
 [207,] 3.252155 3.215063 3.163088 3.070946
 [208,] 3.248832 3.211873 3.160086 3.068276
 [209,] 3.238079 3.201553 3.150372 3.059636
 [210,] 3.197611 3.162714 3.113816 3.027122
 [211,] 3.196471 3.161620 3.112785 3.026206
 [212,] 3.240855 3.204217 3.152880 3.061866
 [213,] 3.231823 3.195548 3.144721 3.054609
 [214,] 3.257129 3.219836 3.167581 3.074942
 [215,] 3.276379 3.238311 3.184971 3.090408
 [216,] 3.261619 3.224146 3.171638 3.078550
 [217,] 3.264186 3.226609 3.173956 3.080612
 [218,] 3.267288 3.229586 3.176758 3.083104
 [219,] 3.275839 3.237793 3.184483 3.089974
 [220,] 3.261036 3.223585 3.171110 3.078080
 [221,] 3.248556 3.211608 3.159836 3.068053
 [222,] 3.258895 3.221531 3.169177 3.076360
 [223,] 3.217515 3.181817 3.131796 3.043114
 [224,] 3.257346 3.220044 3.167777 3.075116
 [225,] 3.232774 3.196462 3.145580 3.055374
 [226,] 3.233536 3.197192 3.146268 3.055986
 [227,] 3.239230 3.202657 3.151412 3.060560
 [228,] 3.236113 3.199666 3.148596 3.058056
 [229,] 3.219017 3.183259 3.133153 3.044321
 [230,] 3.213815 3.178266 3.128453 3.040141
 [231,] 3.172873 3.138972 3.091468 3.007247
 [232,] 3.143054 3.110353 3.064530 2.983288
 [233,] 3.148020 3.115120 3.069017 2.987279
 [234,] 3.125512 3.093518 3.048684 2.969195
 [235,] 3.125541 3.093546 3.048711 2.969218
 [236,] 3.089999 3.059434 3.016603 2.940661
 [237,] 3.117715 3.086034 3.041641 2.962930
 [238,] 3.110692 3.079294 3.035297 2.957288
 [239,] 3.125676 3.093675 3.048832 2.969326
 [240,] 3.129456 3.097302 3.052246 2.972363
 [241,] 3.153600 3.120474 3.074057 2.991761
 [242,] 3.115433 3.083844 3.039579 2.961096
 [243,] 3.102714 3.071637 3.028089 2.950877
 [244,] 3.112779 3.081297 3.037181 2.958964
 [245,] 3.122856 3.090968 3.046284 2.967060
 [246,] 3.122841 3.090954 3.046271 2.967048
 [247,] 3.115490 3.083899 3.039630 2.961142
 [248,] 3.169432 3.135670 3.088360 3.004482
 [249,] 3.151123 3.118097 3.071820 2.989771
 [250,] 3.145183 3.112396 3.066454 2.984999
 [ reached getOption("max.print") -- omitted 831 rows ]

$tau
              [,1] [,2] [,3] [,4]
   [1,] 0.08333333 0.25  0.5    1
   [2,] 0.08333333 0.25  0.5    1
   [3,] 0.08333333 0.25  0.5    1
   [4,] 0.08333333 0.25  0.5    1
   [5,] 0.08333333 0.25  0.5    1
   [6,] 0.08333333 0.25  0.5    1
   [7,] 0.08333333 0.25  0.5    1
   [8,] 0.08333333 0.25  0.5    1
   [9,] 0.08333333 0.25  0.5    1
  [10,] 0.08333333 0.25  0.5    1
  [11,] 0.08333333 0.25  0.5    1
  [12,] 0.08333333 0.25  0.5    1
  [13,] 0.08333333 0.25  0.5    1
  [14,] 0.08333333 0.25  0.5    1
  [15,] 0.08333333 0.25  0.5    1
  [16,] 0.08333333 0.25  0.5    1
  [17,] 0.08333333 0.25  0.5    1
  [18,] 0.08333333 0.25  0.5    1
  [19,] 0.08333333 0.25  0.5    1
  [20,] 0.08333333 0.25  0.5    1
  [21,] 0.08333333 0.25  0.5    1
  [22,] 0.08333333 0.25  0.5    1
  [23,] 0.08333333 0.25  0.5    1
  [24,] 0.08333333 0.25  0.5    1
  [25,] 0.08333333 0.25  0.5    1
  [26,] 0.08333333 0.25  0.5    1
  [27,] 0.08333333 0.25  0.5    1
  [28,] 0.08333333 0.25  0.5    1
  [29,] 0.08333333 0.25  0.5    1
  [30,] 0.08333333 0.25  0.5    1
  [31,] 0.08333333 0.25  0.5    1
  [32,] 0.08333333 0.25  0.5    1
  [33,] 0.08333333 0.25  0.5    1
  [34,] 0.08333333 0.25  0.5    1
  [35,] 0.08333333 0.25  0.5    1
  [36,] 0.08333333 0.25  0.5    1
  [37,] 0.08333333 0.25  0.5    1
  [38,] 0.08333333 0.25  0.5    1
  [39,] 0.08333333 0.25  0.5    1
  [40,] 0.08333333 0.25  0.5    1
  [41,] 0.08333333 0.25  0.5    1
  [42,] 0.08333333 0.25  0.5    1
  [43,] 0.08333333 0.25  0.5    1
  [44,] 0.08333333 0.25  0.5    1
  [45,] 0.08333333 0.25  0.5    1
  [46,] 0.08333333 0.25  0.5    1
  [47,] 0.08333333 0.25  0.5    1
  [48,] 0.08333333 0.25  0.5    1
  [49,] 0.08333333 0.25  0.5    1
  [50,] 0.08333333 0.25  0.5    1
  [51,] 0.08333333 0.25  0.5    1
  [52,] 0.08333333 0.25  0.5    1
  [53,] 0.08333333 0.25  0.5    1
  [54,] 0.08333333 0.25  0.5    1
  [55,] 0.08333333 0.25  0.5    1
  [56,] 0.08333333 0.25  0.5    1
  [57,] 0.08333333 0.25  0.5    1
  [58,] 0.08333333 0.25  0.5    1
  [59,] 0.08333333 0.25  0.5    1
  [60,] 0.08333333 0.25  0.5    1
  [61,] 0.08333333 0.25  0.5    1
  [62,] 0.08333333 0.25  0.5    1
  [63,] 0.08333333 0.25  0.5    1
  [64,] 0.08333333 0.25  0.5    1
  [65,] 0.08333333 0.25  0.5    1
  [66,] 0.08333333 0.25  0.5    1
  [67,] 0.08333333 0.25  0.5    1
  [68,] 0.08333333 0.25  0.5    1
  [69,] 0.08333333 0.25  0.5    1
  [70,] 0.08333333 0.25  0.5    1
  [71,] 0.08333333 0.25  0.5    1
  [72,] 0.08333333 0.25  0.5    1
  [73,] 0.08333333 0.25  0.5    1
  [74,] 0.08333333 0.25  0.5    1
  [75,] 0.08333333 0.25  0.5    1
  [76,] 0.08333333 0.25  0.5    1
  [77,] 0.08333333 0.25  0.5    1
  [78,] 0.08333333 0.25  0.5    1
  [79,] 0.08333333 0.25  0.5    1
  [80,] 0.08333333 0.25  0.5    1
  [81,] 0.08333333 0.25  0.5    1
  [82,] 0.08333333 0.25  0.5    1
  [83,] 0.08333333 0.25  0.5    1
  [84,] 0.08333333 0.25  0.5    1
  [85,] 0.08333333 0.25  0.5    1
  [86,] 0.08333333 0.25  0.5    1
  [87,] 0.08333333 0.25  0.5    1
  [88,] 0.08333333 0.25  0.5    1
  [89,] 0.08333333 0.25  0.5    1
  [90,] 0.08333333 0.25  0.5    1
  [91,] 0.08333333 0.25  0.5    1
  [92,] 0.08333333 0.25  0.5    1
  [93,] 0.08333333 0.25  0.5    1
  [94,] 0.08333333 0.25  0.5    1
  [95,] 0.08333333 0.25  0.5    1
  [96,] 0.08333333 0.25  0.5    1
  [97,] 0.08333333 0.25  0.5    1
  [98,] 0.08333333 0.25  0.5    1
  [99,] 0.08333333 0.25  0.5    1
 [100,] 0.08333333 0.25  0.5    1
 [101,] 0.08333333 0.25  0.5    1
 [102,] 0.08333333 0.25  0.5    1
 [103,] 0.08333333 0.25  0.5    1
 [104,] 0.08333333 0.25  0.5    1
 [105,] 0.08333333 0.25  0.5    1
 [106,] 0.08333333 0.25  0.5    1
 [107,] 0.08333333 0.25  0.5    1
 [108,] 0.08333333 0.25  0.5    1
 [109,] 0.08333333 0.25  0.5    1
 [110,] 0.08333333 0.25  0.5    1
 [111,] 0.08333333 0.25  0.5    1
 [112,] 0.08333333 0.25  0.5    1
 [113,] 0.08333333 0.25  0.5    1
 [114,] 0.08333333 0.25  0.5    1
 [115,] 0.08333333 0.25  0.5    1
 [116,] 0.08333333 0.25  0.5    1
 [117,] 0.08333333 0.25  0.5    1
 [118,] 0.08333333 0.25  0.5    1
 [119,] 0.08333333 0.25  0.5    1
 [120,] 0.08333333 0.25  0.5    1
 [121,] 0.08333333 0.25  0.5    1
 [122,] 0.08333333 0.25  0.5    1
 [123,] 0.08333333 0.25  0.5    1
 [124,] 0.08333333 0.25  0.5    1
 [125,] 0.08333333 0.25  0.5    1
 [126,] 0.08333333 0.25  0.5    1
 [127,] 0.08333333 0.25  0.5    1
 [128,] 0.08333333 0.25  0.5    1
 [129,] 0.08333333 0.25  0.5    1
 [130,] 0.08333333 0.25  0.5    1
 [131,] 0.08333333 0.25  0.5    1
 [132,] 0.08333333 0.25  0.5    1
 [133,] 0.08333333 0.25  0.5    1
 [134,] 0.08333333 0.25  0.5    1
 [135,] 0.08333333 0.25  0.5    1
 [136,] 0.08333333 0.25  0.5    1
 [137,] 0.08333333 0.25  0.5    1
 [138,] 0.08333333 0.25  0.5    1
 [139,] 0.08333333 0.25  0.5    1
 [140,] 0.08333333 0.25  0.5    1
 [141,] 0.08333333 0.25  0.5    1
 [142,] 0.08333333 0.25  0.5    1
 [143,] 0.08333333 0.25  0.5    1
 [144,] 0.08333333 0.25  0.5    1
 [145,] 0.08333333 0.25  0.5    1
 [146,] 0.08333333 0.25  0.5    1
 [147,] 0.08333333 0.25  0.5    1
 [148,] 0.08333333 0.25  0.5    1
 [149,] 0.08333333 0.25  0.5    1
 [150,] 0.08333333 0.25  0.5    1
 [151,] 0.08333333 0.25  0.5    1
 [152,] 0.08333333 0.25  0.5    1
 [153,] 0.08333333 0.25  0.5    1
 [154,] 0.08333333 0.25  0.5    1
 [155,] 0.08333333 0.25  0.5    1
 [156,] 0.08333333 0.25  0.5    1
 [157,] 0.08333333 0.25  0.5    1
 [158,] 0.08333333 0.25  0.5    1
 [159,] 0.08333333 0.25  0.5    1
 [160,] 0.08333333 0.25  0.5    1
 [161,] 0.08333333 0.25  0.5    1
 [162,] 0.08333333 0.25  0.5    1
 [163,] 0.08333333 0.25  0.5    1
 [164,] 0.08333333 0.25  0.5    1
 [165,] 0.08333333 0.25  0.5    1
 [166,] 0.08333333 0.25  0.5    1
 [167,] 0.08333333 0.25  0.5    1
 [168,] 0.08333333 0.25  0.5    1
 [169,] 0.08333333 0.25  0.5    1
 [170,] 0.08333333 0.25  0.5    1
 [171,] 0.08333333 0.25  0.5    1
 [172,] 0.08333333 0.25  0.5    1
 [173,] 0.08333333 0.25  0.5    1
 [174,] 0.08333333 0.25  0.5    1
 [175,] 0.08333333 0.25  0.5    1
 [176,] 0.08333333 0.25  0.5    1
 [177,] 0.08333333 0.25  0.5    1
 [178,] 0.08333333 0.25  0.5    1
 [179,] 0.08333333 0.25  0.5    1
 [180,] 0.08333333 0.25  0.5    1
 [181,] 0.08333333 0.25  0.5    1
 [182,] 0.08333333 0.25  0.5    1
 [183,] 0.08333333 0.25  0.5    1
 [184,] 0.08333333 0.25  0.5    1
 [185,] 0.08333333 0.25  0.5    1
 [186,] 0.08333333 0.25  0.5    1
 [187,] 0.08333333 0.25  0.5    1
 [188,] 0.08333333 0.25  0.5    1
 [189,] 0.08333333 0.25  0.5    1
 [190,] 0.08333333 0.25  0.5    1
 [191,] 0.08333333 0.25  0.5    1
 [192,] 0.08333333 0.25  0.5    1
 [193,] 0.08333333 0.25  0.5    1
 [194,] 0.08333333 0.25  0.5    1
 [195,] 0.08333333 0.25  0.5    1
 [196,] 0.08333333 0.25  0.5    1
 [197,] 0.08333333 0.25  0.5    1
 [198,] 0.08333333 0.25  0.5    1
 [199,] 0.08333333 0.25  0.5    1
 [200,] 0.08333333 0.25  0.5    1
 [201,] 0.08333333 0.25  0.5    1
 [202,] 0.08333333 0.25  0.5    1
 [203,] 0.08333333 0.25  0.5    1
 [204,] 0.08333333 0.25  0.5    1
 [205,] 0.08333333 0.25  0.5    1
 [206,] 0.08333333 0.25  0.5    1
 [207,] 0.08333333 0.25  0.5    1
 [208,] 0.08333333 0.25  0.5    1
 [209,] 0.08333333 0.25  0.5    1
 [210,] 0.08333333 0.25  0.5    1
 [211,] 0.08333333 0.25  0.5    1
 [212,] 0.08333333 0.25  0.5    1
 [213,] 0.08333333 0.25  0.5    1
 [214,] 0.08333333 0.25  0.5    1
 [215,] 0.08333333 0.25  0.5    1
 [216,] 0.08333333 0.25  0.5    1
 [217,] 0.08333333 0.25  0.5    1
 [218,] 0.08333333 0.25  0.5    1
 [219,] 0.08333333 0.25  0.5    1
 [220,] 0.08333333 0.25  0.5    1
 [221,] 0.08333333 0.25  0.5    1
 [222,] 0.08333333 0.25  0.5    1
 [223,] 0.08333333 0.25  0.5    1
 [224,] 0.08333333 0.25  0.5    1
 [225,] 0.08333333 0.25  0.5    1
 [226,] 0.08333333 0.25  0.5    1
 [227,] 0.08333333 0.25  0.5    1
 [228,] 0.08333333 0.25  0.5    1
 [229,] 0.08333333 0.25  0.5    1
 [230,] 0.08333333 0.25  0.5    1
 [231,] 0.08333333 0.25  0.5    1
 [232,] 0.08333333 0.25  0.5    1
 [233,] 0.08333333 0.25  0.5    1
 [234,] 0.08333333 0.25  0.5    1
 [235,] 0.08333333 0.25  0.5    1
 [236,] 0.08333333 0.25  0.5    1
 [237,] 0.08333333 0.25  0.5    1
 [238,] 0.08333333 0.25  0.5    1
 [239,] 0.08333333 0.25  0.5    1
 [240,] 0.08333333 0.25  0.5    1
 [241,] 0.08333333 0.25  0.5    1
 [242,] 0.08333333 0.25  0.5    1
 [243,] 0.08333333 0.25  0.5    1
 [244,] 0.08333333 0.25  0.5    1
 [245,] 0.08333333 0.25  0.5    1
 [246,] 0.08333333 0.25  0.5    1
 [247,] 0.08333333 0.25  0.5    1
 [248,] 0.08333333 0.25  0.5    1
 [249,] 0.08333333 0.25  0.5    1
 [250,] 0.08333333 0.25  0.5    1
 [ reached getOption("max.print") -- omitted 831 rows ]

$r
   [1] 3.550000 3.538842 3.552720 3.541082
   [5] 3.552227 3.537761 3.505029 3.468814
   [9] 3.476318 3.494846 3.462701 3.473242
  [13] 3.507879 3.519130 3.546248 3.520318
  [17] 3.515973 3.533916 3.540816 3.534481
  [21] 3.553977 3.560081 3.544835 3.505558
  [25] 3.515499 3.516216 3.516490 3.504421
  [29] 3.533041 3.527999 3.505725 3.502654
  [33] 3.490419 3.507607 3.494805 3.475603
  [37] 3.484817 3.514903 3.535239 3.560649
  [41] 3.568608 3.577464 3.577094 3.566816
  [45] 3.605194 3.595943 3.613350 3.638628
  [49] 3.624246 3.638777 3.610782 3.595057
  [53] 3.589805 3.598197 3.597013 3.571074
  [57] 3.596525 3.563857 3.562641 3.569732
  [61] 3.580801 3.562728 3.582055 3.588683
  [65] 3.594200 3.604657 3.578938 3.579192
  [69] 3.577289 3.591850 3.578281 3.561481
  [73] 3.580258 3.563779 3.592489 3.592167
  [77] 3.577425 3.608223 3.593507 3.611618
  [81] 3.628826 3.622309 3.648233 3.639412
  [85] 3.657627 3.668200 3.657034 3.624078
  [89] 3.643643 3.673525 3.667862 3.674103
  [93] 3.676690 3.670352 3.640650 3.643894
  [97] 3.654186 3.669084 3.633133 3.610372
 [101] 3.590100 3.595519 3.574625 3.539614
 [105] 3.533815 3.521395 3.527278 3.546542
 [109] 3.570951 3.559714 3.551716 3.549813
 [113] 3.559252 3.587389 3.579009 3.582955
 [117] 3.561019 3.523332 3.529664 3.517851
 [121] 3.490291 3.494454 3.508813 3.484013
 [125] 3.482539 3.493938 3.498490 3.523321
 [129] 3.523696 3.503439 3.487178 3.466753
 [133] 3.455621 3.475253 3.500840 3.497392
 [137] 3.488257 3.514911 3.513224 3.501202
 [141] 3.492953 3.460520 3.453726 3.527524
 [145] 3.534896 3.525337 3.529129 3.485098
 [149] 3.489284 3.478761 3.502800 3.478427
 [153] 3.474560 3.473498 3.463993 3.476991
 [157] 3.489114 3.473059 3.478289 3.482304
 [161] 3.481841 3.456237 3.466647 3.430073
 [165] 3.437459 3.423877 3.398695 3.408769
 [169] 3.360588 3.383209 3.381553 3.387869
 [173] 3.354118 3.369885 3.341591 3.325952
 [177] 3.295729 3.289259 3.313657 3.308392
 [181] 3.309249 3.346025 3.303777 3.284521
 [185] 3.279526 3.285210 3.295000 3.283311
 [189] 3.281714 3.279950 3.255354 3.265252
 [193] 3.278866 3.243953 3.225195 3.241196
 [197] 3.223731 3.213021 3.257231 3.255589
 [201] 3.265181 3.245946 3.248625 3.245135
 [205] 3.270786 3.305003 3.271481 3.268088
 [209] 3.257109 3.215793 3.214628 3.259943
 [213] 3.250722 3.276559 3.296212 3.281144
 [217] 3.283764 3.286931 3.295661 3.280547
 [221] 3.267806 3.278362 3.236114 3.276780
 [225] 3.251693 3.252471 3.258284 3.255102
 [229] 3.237648 3.232336 3.190536 3.160091
 [233] 3.165161 3.142181 3.142211 3.105923
 [237] 3.134220 3.127050 3.142348 3.146207
 [241] 3.170858 3.131890 3.118904 3.129180
 [245] 3.139469 3.139453 3.131948 3.187022
 [249] 3.168329 3.162264 3.152479 3.132748
 [253] 3.119625 3.069987 3.050819 3.061804
 [257] 3.042818 3.018123 3.034240 3.056229
 [261] 3.061187 3.066987 3.066630 3.077889
 [265] 3.061758 3.060031 3.049333 3.067570
 [269] 3.083712 3.076777 3.094727 3.094308
 [273] 3.117604 3.147173 3.139975 3.136624
 [277] 3.127450 3.094094 3.091978 3.109603
 [281] 3.110262 3.106145 3.094882 3.090942
 [285] 3.102621 3.108990 3.102712 3.109262
 [289] 3.101948 3.140572 3.129236 3.109153
 [293] 3.099763 3.115735 3.102565 3.118294
 [297] 3.131082 3.127067 3.122675 3.140737
 [301] 3.150672 3.144269 3.160267 3.136845
 [305] 3.135836 3.172743 3.162823 3.172713
 [309] 3.176647 3.198064 3.193872 3.233750
 [313] 3.228006 3.221467 3.169976 3.153952
 [317] 3.186829 3.222241 3.219652 3.217614
 [321] 3.245700 3.238724 3.220708 3.243536
 [325] 3.202006 3.177719 3.161742 3.142338
 [329] 3.147589 3.145200 3.143190 3.153793
 [333] 3.144260 3.160783 3.144645 3.142982
 [337] 3.120793 3.132457 3.147241 3.101957
 [341] 3.076776 3.094021 3.084982 3.083671
 [345] 3.072273 3.061542 3.064461 3.066496
 [349] 3.076061 3.114669 3.127079 3.154439
 [353] 3.149093 3.142713 3.165199 3.157469
 [357] 3.158190 3.153589 3.139000 3.147841
 [361] 3.140424 3.118107 3.123454 3.117100
 [365] 3.104586 3.103292 3.076296 3.072030
 [369] 3.073758 3.049640 3.044320 3.041834
 [373] 3.028797 3.042626 3.049131 3.035058
 [377] 3.029424 3.017409 3.052834 3.055358
 [381] 3.078135 3.044495 3.039039 3.053225
 [385] 3.066969 3.061752 3.085014 3.071139
 [389] 3.075178 3.069935 3.067051 3.055324
 [393] 3.055936 3.038679 3.044491 3.053471
 [397] 3.080897 3.106037 3.102653 3.111401
 [401] 3.103559 3.082527 3.116471 3.153950
 [405] 3.168777 3.167493 3.160724 3.156570
 [409] 3.146208 3.139304 3.130092 3.133638
 [413] 3.144708 3.147726 3.133236 3.118741
 [417] 3.115866 3.105960 3.102232 3.122851
 [421] 3.099514 3.119025 3.133114 3.152185
 [425] 3.123098 3.110923 3.078176 3.088537
 [429] 3.113651 3.127372 3.124474 3.106189
 [433] 3.112305 3.107770 3.095351 3.075943
 [437] 3.059706 3.049427 3.059503 3.067712
 [441] 3.069884 3.088228 3.064402 3.059869
 [445] 3.048267 3.063199 3.061664 3.054330
 [449] 3.058681 3.058235 3.054319 3.046859
 [453] 3.041180 3.034822 3.061436 3.049763
 [457] 3.081342 3.062073 3.047298 3.057570
 [461] 3.054351 3.049642 3.039422 3.038098
 [465] 3.042358 3.037300 3.037580 3.000909
 [469] 3.014076 3.004481 2.975367 2.980794
 [473] 2.984025 2.999224 2.986465 2.969569
 [477] 3.000562 2.989343 2.988557 2.978036
 [481] 2.980177 2.987740 2.995917 2.976284
 [485] 2.956840 2.955123 2.950071 2.973542
 [489] 2.981874 2.978320 3.006330 3.002090
 [493] 2.975521 2.979392 2.995446 2.996016
 [497] 2.968278 2.928555 2.936433 2.967057
 [501] 2.949975 2.937423 2.909714 2.876826
 [505] 2.855806 2.865701 2.865649 2.839555
 [509] 2.816853 2.829056 2.828561 2.848877
 [513] 2.831983 2.829028 2.821790 2.814391
 [517] 2.804791 2.824961 2.813735 2.823639
 [521] 2.816231 2.783414 2.817983 2.802083
 [525] 2.797659 2.784107 2.817207 2.819790
 [529] 2.791496 2.751494 2.736678 2.773438
 [533] 2.733525 2.739849 2.751518 2.718717
 [537] 2.725465 2.699264 2.703799 2.712065
 [541] 2.719539 2.697862 2.694359 2.739841
 [545] 2.728472 2.729173 2.733591 2.772786
 [549] 2.777655 2.776133 2.799716 2.760730
 [553] 2.754681 2.783530 2.813677 2.824143
 [557] 2.855787 2.870687 2.863166 2.838197
 [561] 2.844385 2.822802 2.817263 2.810491
 [565] 2.764741 2.772507 2.779928 2.775654
 [569] 2.766507 2.739090 2.723725 2.695975
 [573] 2.704232 2.706104 2.694799 2.692381
 [577] 2.677638 2.721344 2.750713 2.729524
 [581] 2.743528 2.713842 2.702755 2.713293
 [585] 2.743462 2.769645 2.766426 2.758835
 [589] 2.762716 2.781753 2.757075 2.762679
 [593] 2.765286 2.737350 2.727367 2.718490
 [597] 2.734534 2.722822 2.730396 2.706467
 [601] 2.710205 2.688809 2.681616 2.678362
 [605] 2.679091 2.663196 2.695600 2.672522
 [609] 2.691634 2.652835 2.681773 2.686395
 [613] 2.711409 2.711479 2.726080 2.702571
 [617] 2.720876 2.746142 2.729759 2.704866
 [621] 2.699013 2.703378 2.722000 2.698232
 [625] 2.712587 2.724420 2.709716 2.704101
 [629] 2.704259 2.683042 2.680379 2.685217
 [633] 2.669877 2.656816 2.666095 2.653054
 [637] 2.647532 2.651215 2.614585 2.603311
 [641] 2.581018 2.607717 2.613501 2.593581
 [645] 2.601158 2.606431 2.594266 2.552843
 [649] 2.581749 2.575985 2.555129 2.561848
 [653] 2.593774 2.616710 2.609371 2.626748
 [657] 2.645852 2.655281 2.634948 2.650107
 [661] 2.670021 2.681284 2.683535 2.681664
 [665] 2.683473 2.677894 2.692658 2.715258
 [669] 2.722244 2.720116 2.727714 2.742871
 [673] 2.752266 2.776025 2.769817 2.755429
 [677] 2.757392 2.782548 2.764268 2.765048
 [681] 2.742421 2.731769 2.716674 2.683446
 [685] 2.692140 2.688139 2.675533 2.669369
 [689] 2.667225 2.674058 2.670222 2.714302
 [693] 2.737054 2.724498 2.736202 2.731355
 [697] 2.682167 2.687308 2.682847 2.685321
 [701] 2.684424 2.673791 2.681222 2.677487
 [705] 2.664810 2.642795 2.648333 2.657305
 [709] 2.656975 2.623279 2.658545 2.688641
 [713] 2.695580 2.698351 2.703730 2.680236
 [717] 2.657811 2.661598 2.639449 2.638964
 [721] 2.649289 2.647645 2.664386 2.643260
 [725] 2.647696 2.605298 2.610075 2.564592
 [729] 2.566204 2.562332 2.539019 2.548684
 [733] 2.551548 2.537893 2.546772 2.545423
 [737] 2.527770 2.533057 2.552726 2.563281
 [741] 2.568925 2.557774 2.527588 2.507565
 [745] 2.526259 2.539673 2.551968 2.534615
 [749] 2.552874 2.549068 2.561530 2.570164
 [753] 2.565552 2.593283 2.567948 2.558190
 [757] 2.554758 2.542257 2.565900 2.566112
 [761] 2.553051 2.555879 2.607319 2.629677
 [765] 2.637671 2.656839 2.665417 2.684225
 [769] 2.686361 2.708768 2.700597 2.719420
 [773] 2.699274 2.695315 2.694450 2.711965
 [777] 2.720171 2.737522 2.741073 2.747178
 [781] 2.748263 2.759628 2.737542 2.762719
 [785] 2.747036 2.711475 2.751927 2.747364
 [789] 2.745357 2.767136 2.760496 2.744809
 [793] 2.751439 2.764981 2.793119 2.796019
 [797] 2.798809 2.818128 2.837095 2.840698
 [801] 2.851617 2.834849 2.830769 2.823625
 [805] 2.798962 2.777674 2.749314 2.742179
 [809] 2.757637 2.754033 2.745782 2.753266
 [813] 2.737429 2.745140 2.723665 2.735331
 [817] 2.737664 2.772816 2.784271 2.779346
 [821] 2.771202 2.777760 2.776227 2.757810
 [825] 2.772499 2.765642 2.764904 2.767490
 [829] 2.785753 2.779041 2.795074 2.824386
 [833] 2.830205 2.859133 2.843249 2.820439
 [837] 2.841955 2.847626 2.840162 2.848081
 [841] 2.845383 2.824779 2.833794 2.819607
 [845] 2.824353 2.819331 2.826142 2.795490
 [849] 2.824270 2.821220 2.805920 2.827396
 [853] 2.838493 2.835296 2.806456 2.795527
 [857] 2.804299 2.796369 2.815403 2.815215
 [861] 2.793929 2.801768 2.797053 2.779592
 [865] 2.767325 2.779812 2.783636 2.787206
 [869] 2.777381 2.796122 2.811228 2.820750
 [873] 2.833513 2.824867 2.807388 2.809464
 [877] 2.827966 2.839922 2.837756 2.835559
 [881] 2.837391 2.831086 2.862453 2.899569
 [885] 2.906674 2.892675 2.870838 2.859669
 [889] 2.867400 2.888025 2.880871 2.879087
 [893] 2.853316 2.824731 2.848023 2.867468
 [897] 2.856393 2.852702 2.840926 2.827492
 [901] 2.803940 2.825222 2.799841 2.809826
 [905] 2.817016 2.830048 2.818453 2.819942
 [909] 2.818768 2.836246 2.812643 2.828938
 [913] 2.857766 2.820510 2.818899 2.844770
 [917] 2.808532 2.821166 2.827941 2.848692
 [921] 2.840973 2.813930 2.781262 2.793097
 [925] 2.814549 2.797640 2.779057 2.776604
 [929] 2.793701 2.815212 2.782690 2.816403
 [933] 2.785483 2.777491 2.777611 2.761836
 [937] 2.753294 2.739333 2.755955 2.734312
 [941] 2.722358 2.719891 2.746614 2.718820
 [945] 2.725283 2.712693 2.691864 2.693476
 [949] 2.690050 2.673233 2.700186 2.704926
 [953] 2.675693 2.681928 2.672122 2.663053
 [957] 2.603618 2.617606 2.619839 2.632373
 [961] 2.631898 2.648560 2.633705 2.633756
 [965] 2.626653 2.634867 2.666574 2.670859
 [969] 2.646292 2.653195 2.652758 2.640057
 [973] 2.634044 2.637998 2.642126 2.671018
 [977] 2.652545 2.681572 2.690554 2.688577
 [981] 2.702903 2.701590 2.722306 2.732686
 [985] 2.721884 2.745666 2.735561 2.733187
 [989] 2.741005 2.734311 2.743407 2.742716
 [993] 2.754991 2.752031 2.754660 2.755123
 [997] 2.746619 2.742952 2.750335 2.719170
 [ reached getOption("max.print") -- omitted 81 entries ]

> plot(bond.vasicek)
Error in x(x) : argument "maturities" is missing, with no default
> library("fExoticOptions", lib.loc="~/R/win-library/3.5")
Warning message:
package ‘fExoticOptions’ was built under R version 3.5.1
> library(fExoticOptions)
> a <- GBSOption("c", 100, 100, 1, 0.02, -0.02, 0.3, title = NULL,
+                description = NULL)
> (z <- a@price)
[1] 10.62678
> a <- GeometricAverageRateOption("c", 100, 100, 1, 0.02, -0.02, 0.3,
+                                 title = NULL, description = NULL)
> (z <- a@price)
[1] 5.889822
> a <- StandardBarrierOption("cuo", 100, 90, 130, 0, 1, 0.02, -0.02, 0.30,
+                            title = NULL, description = NULL)
> x <- a@price
> b <- StandardBarrierOption("cui", 100, 90, 130, 0, 1, 0.02, -0.02, 0.30,
+                            title = NULL, description = NULL)
> y <- b@price
> c <- GBSOption("c", 100, 90, 1, 0.02, -0.02, 0.3, title = NULL,
+                description = NULL)
> z <- c@price
> v <- z - x - y
> v
[1] 0
> vanilla <- GBSOption(TypeFlag = "c", S = 100, X = 90, Time = 1,
+                      r = 0.02, b = -0.02, sigma = 0.3)
> KO <- sapply(100:300, FUN = StandardBarrierOption, TypeFlag = "cuo",
+              S = 100, X = 90, K = 0, Time = 1, r = 0.02, b = -0.02, sigma = 0.30)
Warning messages:
1: In if (X >= H) { :
  the condition has length > 1 and only the first element will be used
2: In if (X < H) { :
  the condition has length > 1 and only the first element will be used
> plot(KO[[1]]@price, type = "l",
+      xlab = "barrier distance from spot",
+      ylab = "price of option",
+      main = "Price of KO converges to plain vanilla")
> abline(h = vanilla@price, col = "red")
> library("plot3D", lib.loc="~/R/win-library/3.5")
Warning message:
package ‘plot3D’ was built under R version 3.5.2
> library(plot3D)
> BS_surface <- function(S, Time, FUN, ...) {
+     require(plot3D)
+     n <- length(S)
+     k <- length(Time)
+     m <- matrix(0, n, k)
+     for (i in 1:n){
+         for (j in 1:k){
+             l <- list(S = S[i], Time = Time[j], ...)
+             m[i,j] <- max(do.call(FUN, l)@price, 0)
+         }
+     }
+ persp3D(z = m, xlab = "underlying", ylab = "Remaining time",
+         zlab = "option price", phi = 30, theta = 20, bty = "b2")
+ }
> BS_surface(seq(1, 200,length = 200), seq(0, 2, length = 200),
+            GBSOption, TypeFlag = "c", X = 90, r = 0.02, b = 0, sigma = 0.3)
> BS_surface(seq(1,200,length = 200), seq(0, 2, length = 200),
+            StandardBarrierOption, TypeFlag = "cuo", H = 130, X = 90, K = 0,
+            r = 0.02, b = -0.02, sigma = 0.30)
> GetGreeks <- function(FUN, arg, epsilon,...) {
+     all_args1 <- all_args2 <- list(...)
+     all_args1[[arg]] <- as.numeric(all_args1[[arg]] + epsilon)
+     all_args2[[arg]] <- as.numeric(all_args2[[arg]] - epsilon)
+     (do.call(FUN, all_args1)@price -
+             do.call(FUN, all_args2)@price) / (2 * epsilon)
+ }
> x <- seq(10, 200, length = 200)
> delta <- vega <- theta <- rho <- rep(0, 200)
> for(i in 1:200){
+     delta[i] <- GetGreeks(FUN = FloatingStrikeLookbackOption,
+                           arg = 2, epsilon = 0.01, "p", x[i], 100, 1, 0.02, -0.02, 0.2)
+     vega[i] <- GetGreeks(FUN = FloatingStrikeLookbackOption,
+                          arg = 7, epsilon = 0.0005, "p", x[i], 100, 1, 0.02, -0.02,
+                          0.2)
+     theta[i] <- GetGreeks(FUN = FloatingStrikeLookbackOption,
+                           arg = 4, epsilon = 1/365, "p", x[i], 100, 1, 0.02, -0.02,
+                           0.2)
+     rho[i] <- GetGreeks(FUN = FloatingStrikeLookbackOption,
+                         arg = 5, epsilon = 0.0001, "p", x[i], 100, 1, 0.02, -0.02, 0.2)
+ }
> par(mfrow = c(2, 2))
> plot(x, delta, type = "l", xlab = "S", ylab = "", main = "Delta")
> plot(x, vega, type = "l", xlab = "S", ylab = "", main = "Vega")
> plot(x, theta, type = "l", xlab = "S", ylab = "", main = "Theta")
> plot(x, rho, type = "l", xlab = "S", ylab = "", main = "Rho")
> dnt1 <- function(S, K, U, L, sigma, T, r, b, N = 20, ploterror = FALSE){
+     if ( L > S | S > U) return(0)
+     Z <- log(U/L)
+     alpha <- -1/2*(2*b/sigma^2 - 1)
+     beta <- -1/4*(2*b/sigma^2 - 1)^2 - 2*r/sigma^2
+     v <- rep(0, N)
+     for (i in 1:N)
+         v[i] <- 2*pi*i*K/(Z^2) * (((S/L)^alpha - (-1)^i*(S/U)^alpha ) /(alpha^2+(i*pi/Z)^2)) * sin(i*pi/Z*log(S/L)) *
+         exp(-1/2 * ((i*pi/Z)^2-beta) * sigma^2*T)
+     if (ploterror) barplot(v, main = "Formula Error");
+     sum(v)
+ }
> print(dnt1(100, 10, 120, 80, 0.1, 0.25, 0.05, 0.03, 20, TRUE))
[1] 9.871619
> print(dnt1(100, 10, 120, 80, 0.03, 0.25, 0.05, 0.03, 50, TRUE))
[1] 9.875778
> dnt1 <- function(S, K, U, L, sigma, Time, r, b) {
+     if ( L > S | S > U) return(0)
+     Z <- log(U/L)
+     alpha <- -1/2*(2*b/sigma^2 - 1)
+     beta <- -1/4*(2*b/sigma^2 - 1)^2 - 2*r/sigma^2
+     p <- 0
+     i <- a <- 1
+     while (abs(a) > 0.0001){
+         a <- 2*pi*i*K/(Z^2) * (((S/L)^alpha - (-1)^i*(S/U)^alpha ) /
+                                    (alpha^2 + (i *pi / Z)^2) ) * sin(i * pi / Z * log(S/L)) *
+             exp(-1/2*((i*pi/Z)^2-beta) * sigma^2 * Time)
+         p <- p + a
+         i <- i + 1
+     }
+     p
+ }
> x <- seq(0.92, 0.96, length = 2000)
> y <- z <- rep(0, 2000)
> for (i in 1:2000){
+     y[i] <- dnt1(x[i], 1e6, 0.96, 0.92, 0.06, 0.25, 0.0025, -0.0250)
+     z[i] <- dnt1(x[i], 1e6, 0.96, 0.92, 0.065, 0.25, 0.0025, -0.0250)
+ }
> matplot(x, cbind(y,z), type = "l", lwd = 2, lty = 1,
+         main = "Price of a DNT with volatility 6% and 6.5%
+         ", cex.main = 0.8, xlab = "Price of underlying" )
> GetGreeks <- function(FUN, arg, epsilon,...) {
+     all_args1 <- all_args2 <- list(...)
+     all_args1[[arg]] <- as.numeric(all_args1[[arg]] + epsilon)
+     all_args2[[arg]] <- as.numeric(all_args2[[arg]] - epsilon)
+     (do.call(FUN, all_args1) -
+             do.call(FUN, all_args2)) / (2 * epsilon)
+ }
> Gamma <- function(FUN, epsilon, S, ...) {
+     arg1 <- list(S, ...)
+     arg2 <- list(S + 2 * epsilon, ...)
+     arg3 <- list(S - 2 * epsilon, ...)
+     y1 <- (do.call(FUN, arg2) - do.call(FUN, arg1)) / (2 * epsilon)
+     y2 <- (do.call(FUN, arg1) - do.call(FUN, arg3)) / (2 * epsilon)
+     (y1 - y2) / (2 * epsilon)
+ }
>
> x = seq(0.9202, 0.9598, length = 200)
> delta <- vega <- theta <- gamma <- rep(0, 200)
> for(i in 1:200){
+     delta[i] <- GetGreeks(FUN = dnt1, arg = 1, epsilon = 0.0001,
+                           x[i], 1000000, 0.96, 0.92, 0.06, 0.5, 0.02, -0.02)
+     vega[i] <- GetGreeks(FUN = dnt1, arg = 5, epsilon = 0.0005,
+                          x[i], 1000000, 0.96, 0.92, 0.06, 0.5, 0.0025, -0.025)
+     theta[i] <- - GetGreeks(FUN = dnt1, arg = 6, epsilon = 1/365,
+                             x[i], 1000000, 0.96, 0.92, 0.06, 0.5, 0.0025, -0.025)
+     gamma[i] <- Gamma(FUN = dnt1, epsilon = 0.0001, S = x[i], K =
+                           1e6, U = 0.96, L = 0.92, sigma = 0.06, Time = 0.5, r = 0.02, b =
+                           -0.02)
+ }
> windows()
> plot(x, vega, type = "l", xlab = "S",ylab = "", main = "Vega")
> ()



Latin Cubes

What is Latin Cubes? The practical applications of Latin crops and related designs are factorial experiment.factorial experiments for me...