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)

finance

option trading


trading

option

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)

Post a Comment

0 Comments