I have recently developed two packages that can accompany the modeling simulation platform I described in my book Business Case Analysis with R: Simulation Tutorials to Support Complex Business Decisions (available at Springer-Nature/Apress and Amazon). These packages are:
For issues, suggestions, or kudos related to the leonRdo and inteRest R packages, please send me a note at [email protected] or visit my Contact web page to send me a note via an easy to use form. I would love to know more about how you might be using these packages.
The leonRdo package reproduces the standard Monte Carlo distributions using the median Latin hypercube sampling (MLHS) process. This process is used to preserve the theoretical shape and mean of a distribution with fewer samples than required using simple Monte Carlo. Most of the distributions are derived from the base R distributions qxxx. Each distribution returns a vector of values of length equal to the number of indicated samples. As a free-be, this package also provides a means to calculate the cross product of two vectors in R3 space and it provides the Iman-Conover correlation method for imposing correlation to a matrix of independent samples via a target correlation matrix.
The full list of functions available in the package is
To install this package, paste the following code in your R console. Don’t forget to specify the ‘lib’ parameter if necessary.
install.packages("https://www.incitedecisiontech.com/packages/leonRdo_0.1.5.tgz",
repos = NULL,
type = "binary")
install.packages("https://www.incitedecisiontech.com/packages/leonRdo_0.1.5.tar.gz",
repos = NULL,
type = "source")
To observe the usefulness of MLHS, consider the distribution of means we obtain by making mutliple independent runs (say, 100) of a Normal distribution (mean = 2, sd = 1) in base R using 1000 trials on each run.
runs <- 100
trials <- 1000
test.means.base <- sapply(1:runs, function(r) mean(rnorm(n = trials, 2, 1)))
Now, observe this distribution of means using MLHS, but this time use just 100 trials.
trials <- 100
test.means.mlhs <- sapply(1:runs, function(r) mean(rnorm_mlhs(n = trials, 2, 1)))
The tradeoff in the much improved stability of the mean of a Normal between trials is the compute-time required to handle the underlying processing of the MLHS. But as we can see, for this example the precision of the base Normal distribution is 684x wider than that of the MLHS while using 10x more trials. To visually reinforce just how stable the MLHS simulation is while using an order of magnitude fewer trials than the base R simulation of the same distribution, consider the scatter plot of the MLHS test means to the base R test means.
The leonRdo package also provides a means to impose correlation on a set of independent trials for a given target correlation matrix. Suppose we have x1
and x2
independent of each other, normally distributed over 500 trials.
x1 <- rnorm_mlhs(500, 0, 1)
x2 <- rnorm_mlhs(500, -1, 2)
Plotting these against each other demonstrates that they are essentially independent.
And using the Pearson correlation we confirm this numerically, too: cor(x1, x2, method = "pearson") =
-0.0223438.
Suppose we want to impose a correlation of 0.75 between these variables. First, set up a matrix to contain the variable values.
V <- matrix(c(x1, x2), ncol = 2)
Then create the desired correlation matrix.
my.corr <- matrix(c(1, 0.75, 0.75, 1), ncol = 2)
## [,1] [,2]
## [1,] 1.00 0.75
## [2,] 0.75 1.00
Of course, this matrix can be imported from a *.csv file or an Excel file.
Now we use the correlate_vars function…
V.corr <- correlate_vars(V, my.corr)
and note that the row values in V have been reordered to obtain the desired approximate correlation.
cor(V.corr, method = "pearson") =
## [,1] [,2]
## [1,] 1.0000000 0.7478373
## [2,] 0.7478373 1.0000000
Since this method uses rank correlation, it does not require that the correlated distributions belong to the same family. Try using another distribution type for one of the independent sample sets while the first remains Normal. Next, try correlating more than two independent sample sets. Remember, though, that your correlation matrix must be positive definite for this method to work.
The inteRest package calculates various time value of money functions such as net present value, internal rate of return, future worth of an annuity, etc.
The full list of functions available in the package is
To install this package, paste the following code in your R console. Don’t forget to specify the ‘lib’ parameter if necessary.
install.packages("https://www.incitedecisiontech.com/packages/inteRest_1.0.tgz",
repos = NULL,
type = "binary")
install.packages("https://www.incitedecisiontech.com/packages/inteRest_1.0.tar.gz",
repos = NULL,
type = "source")
To highlight a few examples of the inteRest package, consider first a cash flow in time.
my.cash.flow <- c(-100, 10, 15, 20, 25, 25, 25, 25, 20, 15)
Assuming an annual discount rate of 0.1, the net present value of the cash flow is found by…
npv <- CalcNPV(my.cash.flow, 0.1)
## [1] 10.67698
You can format the result in dollars…
FormatDollar(npv, cents = TRUE)
## [1] "$10.68"
If an internal rate of return is computable, you can find the IRR by…
irr <- CalcIRR(my.cash.flow, r0 = 0.1, tol = 1e-06, maxiter = 200)
## [1] 0.1263388
Another useful set of calculations are those associated with mortgaging a home. Suppose you want to buy a house for $250,000, financed over 30 years at 4.25% annual interest. Your monthly payment will be…
loan <- 250000
ann.int <- 0.0425
mo.int <- CalcEqPerIR(ann.int, 12)
payment <- CalcPerPayment(L = loan, i = mo.int, N = 360)
## [1] 1218.081
If you want to know the interest, principal, and declining balance on a period- to-period basis…
int <- CalcPerIntPayment(L = loan, i = mo.int, M = 360, N = 1:360)
prin <- CalcPerPrincPayment(L = loan, i = mo.int, M = 360, N = 1:360)
bal <- CalcPerRemBal(L = loan, i = mo.int, M = 360, N = 1:360)
## Period Interest Principal Balance
## 1 1 868.623751 349.4575 249650.543
## 2 2 867.409563 350.6717 249299.871
## 3 3 866.191156 351.8901 248947.981
## 4 4 864.968515 353.1127 248594.868
## 5 5 863.741627 354.3396 248240.528
## 6 6 862.510476 355.5708 247884.958
## 7 7 861.275047 356.8062 247528.152
## 8 8 860.035326 358.0459 247170.106
## 9 9 858.791297 359.2899 246810.816
## 10 10 857.542946 360.5383 246450.277
## 11 351 41.524521 1176.5567 10774.683
## 12 352 37.436581 1180.6447 9594.038
## 13 353 33.334437 1184.7468 8409.291
## 14 354 29.218040 1188.8632 7220.428
## 15 355 25.087341 1192.9939 6027.434
## 16 356 20.942289 1197.1389 4830.295
## 17 357 16.782836 1201.2984 3628.997
## 18 358 12.608931 1205.4723 2423.524
## 19 359 8.420523 1209.6607 1213.864
## 20 360 4.217563 1213.8637 0.000