BSAmericanApproxOption

MonteCarloOption

> mc <- MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,
+                        mcLoops = 50, init = TRUE, innovations.gen = sobolInnovations,
+                        path.gen = wienerPath, payoff.calc = arithmeticAsianPayoff,
+                        antithetic = TRUE, standardization = FALSE, trace = TRUE,
+                        scrambling = 2, seed = 4711)

Add caption

Monte Carlo Simulation Path:


Loop: No Error in MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,  :
  object 'sobolInnovations' not found
> sobolInnovations <- function(mcSteps, pathLength, init, ...) {
+     # Create and return Normal Sobol Innovations:
+     rnorm.sobol(mcSteps, pathLength, init, ...)
+ }
> wienerPath <- function(eps) {
+     # Note, the option parameters must be globally defined!
+     # Generate and return the Paths:
+     (b-sigma*sigma/2)*delta.t + sigma*sqrt(delta.t)*eps
+ }
> plainVanillaPayoff <- function(path) {
+     # Note, the option parameters must be globally defined!
+     # Compute the Call/Put Payoff Value:
+     ST <- S*exp(sum(path))
+     if (TypeFlag == "c") payoff <- exp(-r*Time)*max(ST-X, 0)
+     if (TypeFlag == "p") payoff <- exp(-r*Time)*max(0, X-ST)
+     # Return Value:
+     payoff
+ }
> arithmeticAsianPayoff <- function(path) {
+     # Note, the option parameters must be globally defined!
+     # Compute the Call/Put Payoff Value:
+     SM <- mean(S*exp(cumsum(path)))
+     if (TypeFlag == "c") payoff <- exp(-r*Time)*max(SM-X, 0)
+     if (TypeFlag == "p") payoff <- exp(-r*Time)*max(0, X-SM)
+     # Return Value:
+     payoff
+ }
> TypeFlag <- "c"; S <- 100; X <- 100
> Time <- 1/12; sigma <- 0.4; r <- 0.10; b <- 0.1
>
> mc <- MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,
+                        mcLoops = 50, init = TRUE, innovations.gen = sobolInnovations,
+                        path.gen = wienerPath, payoff.calc = arithmeticAsianPayoff,
+                        antithetic = TRUE, standardization = FALSE, trace = TRUE,
+                        scrambling = 2, seed = 4711)

Monte Carlo Simulation Path:


Loop: No
Loop: 1 : 3.253438 3.253438
Loop: 2 : 2.635154 2.944296
Loop: 3 : 2.6733 2.853964
Loop: 4 : 2.748978 2.827718
Loop: 5 : 2.795612 2.821296
Loop: 6 : 2.904628 2.835185
Loop: 7 : 2.76193 2.82472
Loop: 8 : 2.943871 2.839614
Loop: 9 : 3.628506 2.927269
Loop: 10 : 3.549677 2.989509
Loop: 11 : 2.623178 2.956207
Loop: 12 : 3.26091 2.981599
Loop: 13 : 3.206552 2.998903
Loop: 14 : 3.188396 3.012438
Loop: 15 : 3.080368 3.016967
Loop: 16 : 2.469462 2.982748
Loop: 17 : 2.573869 2.958696
Loop: 18 : 2.736678 2.946362
Loop: 19 : 3.275842 2.963703
Loop: 20 : 2.357823 2.933409
Loop: 21 : 2.6898 2.921808
Loop: 22 : 2.780071 2.915366
Loop: 23 : 2.767824 2.908951
Loop: 24 : 3.090034 2.916496
Loop: 25 : 2.520999 2.900676
Loop: 26 : 2.854595 2.898904
Loop: 27 : 3.274827 2.912827
Loop: 28 : 3.588779 2.936968
Loop: 29 : 2.862528 2.934401
Loop: 30 : 2.987253 2.936163
Loop: 31 : 2.726644 2.929404
Loop: 32 : 2.345068 2.911144
Loop: 33 : 2.849574 2.909278
Loop: 34 : 3.391248 2.923453
Loop: 35 : 2.738561 2.918171
Loop: 36 : 3.244874 2.927246
Loop: 37 : 3.337476 2.938333
Loop: 38 : 2.790049 2.934431
Loop: 39 : 2.308579 2.918383
Loop: 40 : 3.427203 2.931104
Loop: 41 : 2.497448 2.920527
Loop: 42 : 2.562182 2.911995
Loop: 43 : 2.821961 2.909901
Loop: 44 : 2.496715 2.900511
Loop: 45 : 2.627413 2.894442
Loop: 46 : 2.676602 2.889706
Loop: 47 : 2.938307 2.89074
Loop: 48 : 3.024983 2.893537
Loop: 49 : 2.443328 2.884349
Loop: 50 : 3.336588 2.893394
> par(mfrow = c(1, 1))
> mcPrice <- cumsum(mc)/(1:length(mc))
> plot(mcPrice, type = "l", main = "Arithmetic Asian Option",
+      xlab = "Monte Carlo Loops", ylab = "Option Price")
>
> if(FALSE) { #   ... requires(fExoticOptions)
+     TW <- TurnbullWakemanAsianApproxOption(
+         TypeFlag = "c", S = 100, SA = 100, X = 100,
+         Time = 1/12, time = 1/12, tau = 0 , r = 0.1,
+         b = 0.1, sigma = 0.4)$price
+     print(TW)
+ } else
+     TW <- 2.859122
>
> abline(h = TW, col = 2)
> RollGeskeWhaleyOption(S = 80, X = 82, time1 = 1/4,
+                       Time2 = 1/3, r = 0.06, D = 4, sigma = 0.30)

Title:
 Roll Geske Whaley Option

Call:
 RollGeskeWhaleyOption(S = 80, X = 82, time1 = 1/4, Time2 = 1/3,
     r = 0.06, D = 4, sigma = 0.3)

Parameters:
       Value:         
 S     80             
 X     82             
 time1 0.25           
 Time2 0.333333333333333
 r     0.06           
 D     4               
 sigma 0.3             

Option Price:
 4.38603

Description:
 Sun Dec 01 22:43:35 2019

> BAWAmericanApproxOption(TypeFlag = "p", S = 100,
+                         X = 100, Time = 0.5, r = 0.10, b = 0, sigma = 0.25)

Title:
 BAW American Approximated Option

Call:
 BAWAmericanApproxOption(TypeFlag = "p", S = 100, X = 100, Time = 0.5,
     r = 0.1, b = 0, sigma = 0.25)

Parameters:
          Value:
 TypeFlag p   
 S        100 
 X        100 
 Time     0.5 
 r        0.1 
 b        0   
 sigma    0.25 

Option Price:
 6.801362

Description:
 Sun Dec 01 22:44:09 2019

> Black76Option(TypeFlag = "c", FT = 100, X = 100,
+               Time = 0.5, r = 0.10, sigma = 0.25)

Title:
 Black 76 Option Valuation

Call:
 Black76Option(TypeFlag = "c", FT = 100, X = 100, Time = 0.5,
     r = 0.1, sigma = 0.25)

Parameters:
          Value:
 TypeFlag c   
 FT       100 
 X        100 
 Time     0.5 
 r        0.1 
 sigma    0.25 

Option Price:
 6.699709

Description:
 Sun Dec 01 22:44:09 2019

> BSAmericanApproxOption(TypeFlag = "c", S = 42, X = 40,
+                        Time = 0.75, r = 0.04, b = 0.04-0.08, sigma = 0.35)

Title:
 BS American Approximated Option

Call:
 BSAmericanApproxOption(TypeFlag = "c", S = 42, X = 40, Time = 0.75,
     r = 0.04, b = 0.04 - 0.08, sigma = 0.35)

Parameters:
             Value:         
 TypeFlag    c             
 S           42             
 X           40             
 Time        0.75           
 r           0.04           
 b           -0.04         
 sigma       0.35           
 TrigerPrice 57.5994499306841

Option Price:
 5.270405

Description:
 Sun Dec 01 22:44:33 2019
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo
Version 0.4-0 included new data defaults. See ?getSymbols.
> library("SkewHyperbolic", lib.loc="~/R/win-library/3.6")
> mc <- MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,
+                        mcLoops = 50, init = TRUE, innovations.gen = sobolInnovations,
+                        path.gen = wienerPath, payoff.calc = arithmeticAsianPayoff,
+                        antithetic = TRUE, standardization = FALSE, trace = TRUE,
+                        scrambling = 2, seed = 4711)

Monte Carlo Simulation Path:


Loop: No Error in MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,  :
  object 'sobolInnovations' not found
> sobolInnovations <- function(mcSteps, pathLength, init, ...) {
+     # Create and return Normal Sobol Innovations:
+     rnorm.sobol(mcSteps, pathLength, init, ...)
+ }
> wienerPath <- function(eps) {
+     # Note, the option parameters must be globally defined!
+     # Generate and return the Paths:
+     (b-sigma*sigma/2)*delta.t + sigma*sqrt(delta.t)*eps
+ }
> plainVanillaPayoff <- function(path) {
+     # Note, the option parameters must be globally defined!
+     # Compute the Call/Put Payoff Value:
+     ST <- S*exp(sum(path))
+     if (TypeFlag == "c") payoff <- exp(-r*Time)*max(ST-X, 0)
+     if (TypeFlag == "p") payoff <- exp(-r*Time)*max(0, X-ST)
+     # Return Value:
+     payoff
+ }
> arithmeticAsianPayoff <- function(path) {
+     # Note, the option parameters must be globally defined!
+     # Compute the Call/Put Payoff Value:
+     SM <- mean(S*exp(cumsum(path)))
+     if (TypeFlag == "c") payoff <- exp(-r*Time)*max(SM-X, 0)
+     if (TypeFlag == "p") payoff <- exp(-r*Time)*max(0, X-SM)
+     # Return Value:
+     payoff
+ }
> TypeFlag <- "c"; S <- 100; X <- 100
> Time <- 1/12; sigma <- 0.4; r <- 0.10; b <- 0.1
>
> mc <- MonteCarloOption(delta.t = 1/360, pathLength = 30, mcSteps = 5000,
+                        mcLoops = 50, init = TRUE, innovations.gen = sobolInnovations,
+                        path.gen = wienerPath, payoff.calc = arithmeticAsianPayoff,
+                        antithetic = TRUE, standardization = FALSE, trace = TRUE,
+                        scrambling = 2, seed = 4711)

Monte Carlo Simulation Path:


Loop: No
Loop: 1 : 3.253438 3.253438
Loop: 2 : 2.635154 2.944296
Loop: 3 : 2.6733 2.853964
Loop: 4 : 2.748978 2.827718
Loop: 5 : 2.795612 2.821296
Loop: 6 : 2.904628 2.835185
Loop: 7 : 2.76193 2.82472
Loop: 8 : 2.943871 2.839614
Loop: 9 : 3.628506 2.927269
Loop: 10 : 3.549677 2.989509
Loop: 11 : 2.623178 2.956207
Loop: 12 : 3.26091 2.981599
Loop: 13 : 3.206552 2.998903
Loop: 14 : 3.188396 3.012438
Loop: 15 : 3.080368 3.016967
Loop: 16 : 2.469462 2.982748
Loop: 17 : 2.573869 2.958696
Loop: 18 : 2.736678 2.946362
Loop: 19 : 3.275842 2.963703
Loop: 20 : 2.357823 2.933409
Loop: 21 : 2.6898 2.921808
Loop: 22 : 2.780071 2.915366
Loop: 23 : 2.767824 2.908951
Loop: 24 : 3.090034 2.916496
Loop: 25 : 2.520999 2.900676
Loop: 26 : 2.854595 2.898904
Loop: 27 : 3.274827 2.912827
Loop: 28 : 3.588779 2.936968
Loop: 29 : 2.862528 2.934401
Loop: 30 : 2.987253 2.936163
Loop: 31 : 2.726644 2.929404
Loop: 32 : 2.345068 2.911144
Loop: 33 : 2.849574 2.909278
Loop: 34 : 3.391248 2.923453
Loop: 35 : 2.738561 2.918171
Loop: 36 : 3.244874 2.927246
Loop: 37 : 3.337476 2.938333
Loop: 38 : 2.790049 2.934431
Loop: 39 : 2.308579 2.918383
Loop: 40 : 3.427203 2.931104
Loop: 41 : 2.497448 2.920527
Loop: 42 : 2.562182 2.911995
Loop: 43 : 2.821961 2.909901
Loop: 44 : 2.496715 2.900511
Loop: 45 : 2.627413 2.894442
Loop: 46 : 2.676602 2.889706
Loop: 47 : 2.938307 2.89074
Loop: 48 : 3.024983 2.893537
Loop: 49 : 2.443328 2.884349
Loop: 50 : 3.336588 2.893394
> par(mfrow = c(1, 1))
> mcPrice <- cumsum(mc)/(1:length(mc))
> plot(mcPrice, type = "l", main = "Arithmetic Asian Option",
+      xlab = "Monte Carlo Loops", ylab = "Option Price")
>
> if(FALSE) { #   ... requires(fExoticOptions)
+     TW <- TurnbullWakemanAsianApproxOption(
+         TypeFlag = "c", S = 100, SA = 100, X = 100,
+         Time = 1/12, time = 1/12, tau = 0 , r = 0.1,
+         b = 0.1, sigma = 0.4)$price
+     print(TW)
+ } else
+     TW <- 2.859122
>
> abline(h = TW, col = 2)
> RollGeskeWhaleyOption(S = 80, X = 82, time1 = 1/4,
+                       Time2 = 1/3, r = 0.06, D = 4, sigma = 0.30)

Title:
 Roll Geske Whaley Option

Call:
 RollGeskeWhaleyOption(S = 80, X = 82, time1 = 1/4, Time2 = 1/3,
     r = 0.06, D = 4, sigma = 0.3)

Parameters:
       Value:         
 S     80             
 X     82             
 time1 0.25           
 Time2 0.333333333333333
 r     0.06           
 D     4               
 sigma 0.3             

Option Price:
 4.38603

Description:
 Sun Dec 01 22:43:35 2019

> BAWAmericanApproxOption(TypeFlag = "p", S = 100,
+                         X = 100, Time = 0.5, r = 0.10, b = 0, sigma = 0.25)

Title:
 BAW American Approximated Option

Call:
 BAWAmericanApproxOption(TypeFlag = "p", S = 100, X = 100, Time = 0.5,
     r = 0.1, b = 0, sigma = 0.25)

Parameters:
          Value:
 TypeFlag p   
 S        100 
 X        100 
 Time     0.5 
 r        0.1 
 b        0   
 sigma    0.25 

Option Price:
 6.801362

Description:
 Sun Dec 01 22:44:09 2019

> Black76Option(TypeFlag = "c", FT = 100, X = 100,
+               Time = 0.5, r = 0.10, sigma = 0.25)

Title:
 Black 76 Option Valuation

Call:
 Black76Option(TypeFlag = "c", FT = 100, X = 100, Time = 0.5,
     r = 0.1, sigma = 0.25)

Parameters:
          Value:
 TypeFlag c   
 FT       100 
 X        100 
 Time     0.5 
 r        0.1 
 sigma    0.25 

Option Price:
 6.699709

Description:
 Sun Dec 01 22:44:09 2019

> BSAmericanApproxOption(TypeFlag = "c", S = 42, X = 40,
+                        Time = 0.75, r = 0.04, b = 0.04-0.08, sigma = 0.35)

Title:
 BS American Approximated Option

Call:
 BSAmericanApproxOption(TypeFlag = "c", S = 42, X = 40, Time = 0.75,
     r = 0.04, b = 0.04 - 0.08, sigma = 0.35)

Parameters:
             Value:         
 TypeFlag    c             
 S           42             
 X           40             
 Time        0.75           
 r           0.04           
 b           -0.04         
 sigma       0.35           
 TrigerPrice 57.5994499306841

Option Price:
 5.270405

Description:
 Sun Dec 01 22:44:33 2019

> CRRTree = BinomialTreeOption(TypeFlag = "pa", S = 50, X = 50,
+                              Time = 0.4167, r = 0.1, b = 0.1, sigma = 0.4, n = 5)
> BinomialTreePlot(CRRTree, dy = 1, cex = 0.8, ylim = c(-6, 7),
+                  xlab = "n", ylab = "Option Value")
> title(main = "Option Tree")
> model = list(lambda = 4, omega = 8e-5, alpha = 6e-5,
+              beta = 0.7, gamma = 0, rf = 0)
> ts = hngarchSim(model = model, n = 500, n.start = 100)
> par(mfrow = c(2, 1), cex = 0.75)
> ts.plot(ts, col = "steelblue", main = "HN Garch Symmetric Model")
> grid()
> mle = hngarchFit(model = model, x = ts, symmetric = TRUE)
Warning message:
In nlm(.llhHNGarch, par.start, trace = trace, symmetric = symmetric,  :
  NA/Inf replaced by maximum positive value
> mle

Title:
 Heston-Nandi Garch Parameter Estimation

Call:
 hngarchFit(x = ts, model = model, symmetric = TRUE)

Parameters:
lambda   omega   alpha    beta   gamma      rf 
 4e+00   8e-05   6e-05   7e-01   0e+00   0e+00 

Coefficients: lambda, omega, alpha, beta, gamma
   lambda      omega      alpha       beta      gamma 
3.586e+00  3.942e-04  4.147e-05  3.271e-07  0.000e+00 

Log-Likelihood:
 2480.428

Persistence and Variance:
 3.270599e-07
 0.0004356557

Description:
 Sun Dec 01 23:48:32 2019 by user: ADMIN

> model = list(lambda = 4, omega = 8e-5, alpha = 6e-5,
+              beta = 0.7, gamma = 0, rf = 0)
> ts = hngarchSim(model = model, n = 500, n.start = 100)
> par(mfrow = c(2, 1), cex = 0.75)
> ts.plot(ts, col = "steelblue", main = "HN Garch Symmetric Model")
> grid()
>
> mle = hngarchFit(model = model, x = ts, symmetric = TRUE)
Warning message:
In nlm(.llhHNGarch, par.start, trace = trace, symmetric = symmetric,  :
  NA/Inf replaced by maximum positive value
> mle

Title:
 Heston-Nandi Garch Parameter Estimation

Call:
 hngarchFit(x = ts, model = model, symmetric = TRUE)

Parameters:
lambda   omega   alpha    beta   gamma      rf 
 4e+00   8e-05   6e-05   7e-01   0e+00   0e+00 

Coefficients: lambda, omega, alpha, beta, gamma
   lambda      omega      alpha       beta      gamma 
2.267e+00  2.270e-04  6.071e-05  4.850e-01  0.000e+00 

Log-Likelihood:
 2419.411

Persistence and Variance:
 0.4849522
 0.0005586993

Description:
 Sun Dec 01 23:49:12 2019 by user: ADMIN

> par(mfrow = c(3, 1), cex = 0.75)
> summary(mle)

Title:
Heston-Nandi Garch Parameter Estimation

Call:
hngarchFit(x = ts, model = model, symmetric = TRUE)

Parameters:
lambda   omega   alpha    beta   gamma      rf 
 4e+00   8e-05   6e-05   7e-01   0e+00   0e+00 

Coefficients: lambda, omega, alpha, beta, gamma
   lambda      omega      alpha       beta      gamma 
2.267e+00  2.270e-04  6.071e-05  4.850e-01  0.000e+00 

Log-Likelihood:
2419.411

Persistence and Variance:
0.4849522
0.0005586993
> hngarchStats(mle$model)
$mean
[1] 0.001266826

$variance
[1] 0.0005587489

$skewness
[1] 0.004965422

$kurtosis
[1] 3.092809

$persistence
[1] 0.4849522

$leverage
[1] 0

$meansigma2
[1] 0.0005586993

$meansigma4
[1] 3.217825e-07

$meansigma6
[1] 1.925689e-10

$meansigma8
[1] 1.209681e-13

> model = list(lambda = -0.5, omega = 2.3e-6, alpha = 2.9e-6,
+              beta = 0.85, gamma = 184.25)
> S = X = 100
> Time.inDays = 252
> r.daily = 0.05/Time.inDays
> sigma.daily = sqrt((model$omega + model$alpha) /
+                        (1 - model$beta - model$alpha * model$gamma^2))
> data.frame(S, X, r.daily, sigma.daily)
    S   X      r.daily sigma.daily
1 100 100 0.0001984127  0.01004349
> HNG = GBS = Diff = NULL
> for (TypeFlag in c("c", "p")) {
+     HNG = c(HNG, HNGOption(TypeFlag, model = model, S = S, X = X,
+                            Time.inDays = Time.inDays, r.daily = r.daily)$price )
+     GBS = c(GBS, GBSOption(TypeFlag, S = S, X = X, Time = Time.inDays,
+                            r = r.daily, b = r.daily, sigma = sigma.daily)@price) }
> Options = cbind(HNG, GBS, Diff = round(100*(HNG-GBS)/GBS, digits=2))
> row.names(Options) <- c("Call", "Put")
> data.frame(Options)
          HNG      GBS Diff
Call 8.992100 8.939049 0.59
Put  4.115042 4.061991 1.31
> HNG = GBS = Diff = NULL
> for (TypeFlag in c("c", "p")) {
+     HNG = c(HNG, HNGOption(TypeFlag, model = model, S = S, X = X,
+                            Time.inDays = Time.inDays, r.daily = r.daily)$price )
+     GBS = c(GBS, GBSOption(TypeFlag, S = S, X = X, Time = Time.inDays,
+                            r = r.daily, b = r.daily, sigma = sigma.daily)@price) }
> Options = cbind(HNG, GBS, Diff = round(100*(HNG-GBS)/GBS, digits=2))
> row.names(Options) <- c("Call", "Put")
> data.frame(Options)
          HNG      GBS Diff
Call 8.992100 8.939049 0.59
Put  4.115042 4.061991 1.31
> Selection = c("Delta", "Gamma")
> HNG = GBS = NULL
> for (i in 1:2){
+     HNG = c(HNG, HNGGreeks(Selection[i], TypeFlag = "c", model = model,
+                            S = 100, X = 100, Time = Time.inDays, r = r.daily) )
+     GBS = c(GBS, GBSGreeks(Selection[i], TypeFlag = "c", S = 100, X = 100,
+                            Time = Time.inDays, r = r.daily, b = r.daily, sigma = sigma.daily) ) }
> Greeks = cbind(HNG, GBS, Diff = round(100*(HNG-GBS)/GBS, digits = 2))
> row.names(Greeks) <- Selection
> data.frame(Greeks)
             HNG        GBS  Diff
Delta 0.67395339 0.65295993  3.22
Gamma 0.02211149 0.02315963 -4.53
> par(mfrow = c(2, 2), cex = 0.75)
> runif.halton(n = 10, dimension = 5)
        [,1]       [,2] [,3]       [,4]       [,5]
 [1,] 0.5000 0.33333333 0.20 0.14285714 0.09090909
 [2,] 0.2500 0.66666667 0.40 0.28571429 0.18181818
 [3,] 0.7500 0.11111111 0.60 0.42857143 0.27272727
 [4,] 0.1250 0.44444444 0.80 0.57142857 0.36363636
 [5,] 0.6250 0.77777778 0.04 0.71428571 0.45454545
 [6,] 0.3750 0.22222222 0.24 0.85714286 0.54545455
 [7,] 0.8750 0.55555556 0.44 0.02040816 0.63636364
 [8,] 0.0625 0.88888889 0.64 0.16326531 0.72727273
 [9,] 0.5625 0.03703704 0.84 0.30612245 0.81818182
[10,] 0.3125 0.37037037 0.08 0.44897959 0.90909091
> hist(runif.halton(n = 5000, dimension = 1), main = "Uniform Halton",
+      xlab = "x", col = "steelblue3", border = "white")
> rnorm.halton(n = 10, dimension = 5)
            [,1]       [,2]       [,3]       [,4]
 [1,]  0.0000000 -0.4307273 -0.8416212 -1.0675705
 [2,] -0.6744898  0.4307273 -0.2533471 -0.5659489
 [3,]  0.6744898 -1.2206404  0.2533471 -0.1800124
 [4,] -1.1503494 -0.1397103  0.8416212  0.1800124
 [5,]  0.3186394  0.7647097 -1.7506861  0.5659489
 [6,] -0.3186394 -0.7647097 -0.7063026  1.0675705
 [7,]  1.1503494  0.1397103 -0.1509692 -2.0453910
 [8,] -1.5341206  1.2206404  0.3584588 -0.9811260
 [9,]  0.1573107 -1.7861556  0.9944579 -0.5068717
[10,] -0.4887764 -0.3308726 -1.4050716 -0.1282398
            [,5]
 [1,] -1.3351778
 [2,] -0.9084579
 [3,] -0.6045854
 [4,] -0.3487557
 [5,] -0.1141853
 [6,]  0.1141853
 [7,]  0.3487557
 [8,]  0.6045854
 [9,]  0.9084579
[10,]  1.3351778
> hist(rnorm.halton(n = 5000, dimension = 1), main = "Normal Halton",
+      xlab = "x", col = "steelblue3", border = "white")
> runif.sobol(n = 10, dimension = 5, scrambling = 3)
             [,1]       [,2]      [,3]      [,4]       [,5]
 [1,] 0.950436354 0.68109632 0.7964227 0.7789565 0.23022611
 [2,] 0.505117714 0.10001624 0.5394356 0.1626085 0.29512087
 [3,] 0.174511582 0.78909695 0.1558686 0.5312862 0.65873003
 [4,] 0.677048266 0.93456304 0.3741792 0.9448425 0.60587275
 [5,] 0.002695477 0.23692979 0.9461251 0.3135361 0.46957344
 [6,] 0.447998971 0.54445767 0.7023838 0.7443066 0.03626848
 [7,] 0.872369766 0.34755683 0.1179214 0.1217073 0.91564041
 [8,] 0.418673068 0.05109646 0.2309225 0.9069471 0.13713892
 [9,] 0.776573598 0.86932755 0.5893980 0.2873443 0.75079370
[10,] 0.706255794 0.41129673 0.8402236 0.6534177 0.69252151
> hist(runif.sobol(5000, 1, scrambling = 2), main = "Uniform Sobol",
+      xlab = "x", col = "steelblue3", border = "white")
> rnorm.sobol(n = 10, dimension = 5, scrambling = 3)
             [,1]       [,2]       [,3]        [,4]
 [1,]  1.64909933  0.4707667  0.8289112  0.76867389
 [2,]  0.01282859 -1.2814591  0.0990120 -0.98379378
 [3,] -0.93648572  0.8032918 -1.0115836  0.07850354
 [4,]  0.45946060  1.5106647 -0.3208048  1.59677855
 [5,] -2.78269455 -0.7162134  1.6083904 -0.48585188
 [6,] -0.13071858  0.1116705  0.5312691  0.65667975
 [7,]  1.13766483 -0.3919250 -1.1854420 -1.16649456
 [8,] -0.20528928 -1.6343141 -0.7358122  1.32218750
 [9,]  0.76067235  1.1232180  0.2259968 -0.56115990
[10,]  0.54247926 -0.2242106  0.9953773  0.39456417
             [,5]
 [1,] -0.73810241
 [2,] -0.53848580
 [3,]  0.40899966
 [4,]  0.26857795
 [5,] -0.07634218
 [6,] -1.79573322
 [7,]  1.37633104
 [8,] -1.09326416
 [9,]  0.67698953
[10,]  0.50301041
> hist(rnorm.sobol(5000, 1, scrambling = 2), main = "Normal Sobol",
+      xlab = "x", col = "steelblue3", border = "white")
>
> runif.pseudo(n = 10, dimension = 5)
            [,1]       [,2]       [,3]       [,4]
 [1,] 0.79259689 0.70324016 0.01059594 0.65264777
 [2,] 0.08176445 0.11419372 0.63854853 0.14695273
 [3,] 0.43849399 0.82706865 0.66147000 0.86251537
 [4,] 0.03803732 0.62726556 0.91649996 0.01114529
 [5,] 0.54376884 0.14618442 0.76001064 0.15951197
 [6,] 0.51031435 0.05059134 0.75049796 0.84453269
 [7,] 0.59332813 0.59546093 0.02149888 0.90612719
 [8,] 0.47104987 0.25792078 0.42134953 0.60137007
 [9,] 0.67590430 0.55467253 0.14884958 0.39733037
[10,] 0.43858293 0.97749860 0.43211798 0.05088842
            [,5]
 [1,] 0.76051934
 [2,] 0.08410827
 [3,] 0.15209821
 [4,] 0.98667148
 [5,] 0.68380329
 [6,] 0.88421234
 [7,] 0.85896586
 [8,] 0.88135507
 [9,] 0.53711930
[10,] 0.49490493
> rnorm.pseudo(n = 10, dimension = 5)
             [,1]       [,2]       [,3]        [,4]
 [1,] -0.04538924  0.6446193 -0.3001541 -0.23976680
 [2,]  0.20136859 -0.3722682  1.4377910  0.92186264
 [3,] -0.40949839 -0.2485865 -1.7438494 -0.07028057
 [4,]  1.19562289 -1.1394959  1.0015495 -0.31014755
 [5,] -1.20001466 -0.1843586 -1.1581102  0.19796954
 [6,]  1.60942857  0.4281002  1.9370910  0.15939101
 [7,]  2.03454796  1.2277465  0.7349689 -0.20734738
 [8,] -1.03831513  1.0760711 -0.6159524 -1.25480079
 [9,] -0.07567954  0.8022126  0.1477066  1.56488521
[10,] -0.95910380 -0.8359269 -1.2470675 -0.19151533
            [,5]
 [1,]  0.6009371
 [2,]  0.5307484
 [3,] -1.1003838
 [4,] -0.8175360
 [5,] -0.1659358
 [6,] -0.9424639
 [7,]  0.6323741
 [8,] -1.6318668
 [9,]  1.1514525
[10,] -1.4720757