ST303 Stochastic Simumlation

Class 2

Question 1

The inverse tangent function arctan is usually calculated by \[ \arctan(x) = \int_0^x \frac{1}{1 + y^2} dy \,. \] Write a function that uses Monte Carlo simulation to approximate arctan. You can use the same sample for all values of \(x\) or different samples; the latter is much slower.

Use the function to plot the function for values between \(1\) and \(2\) (calculate it for \(1.01, 1.02, 1.03\) and so on till \(2\)).

In particular, we are interested in \(x = \sqrt{3}\) (we know that this is \(\pi/ 3\)).

Compare your results with R function atan.

Question 1

Transformation: \(z = \frac{y}{x}\)

• \(y = z x\)
• \(y'(z) = x\)
• \(z(0) = 0\), \(z(x) = 1\)

\[ \begin{aligned} \arctan(x) &= \int_0^x \frac{1}{1 + y^2} dy \\ &= \int_0^1 \frac{x}{1 + (zx)^2} dz \\ &= \textrm{E}_{\mathcal{U}(0,1)}\left[\frac{x}{1 + (Ux)^2}\right] \\ &\approx \frac{1}{n} \sum_{i = 1}^n \frac{x}{1 + (U_i x)^2}, \quad U_i \overset{\textrm{i.i.d.}}{\sim} \mathcal{U}(0,1) \end{aligned} \]

Question 1

\[ \arctan(x) = \int_0^x \frac{1}{1 + y^2} dy \,. \]

Shortcut

We can think of \(y\) as \(Y \sim \mathcal{U}(0,x)\) with pdf \(f(x) = 1 / x\).

Then

\[ \begin{aligned} \arctan(x) &= \int_0^x \frac{1}{1 + y^2} dy \\ &= \int_0^1 \frac{x}{1 + y^2} \frac{1}{x} dz \\ &= \textrm{E}_{\mathcal{U}(0,x)}\left[\frac{x}{1 + U^2}\right] \\ &\approx \frac{1}{n} \sum_{i = 1}^n \frac{x}{1 + (U_i x)^2}, \quad U_i \overset{\textrm{i.i.d.}}{\sim} \mathcal{U}(0,1) \end{aligned} \]

Question 1

atan_approx <- function(x, u = NULL, seed = 123, N = 10^5){
    if(is.null(u)){
        set.seed(seed)
        u <- runif(N) 
    }
    z <- x/(1+u^2*x^2)
    return (mean(z))
}

set.seed(123)
N <- 10^5
u <- runif(N)
x <- sqrt(3)

cat("MC approx. of arctan(x^(0.5)): ", atan_approx(x, u), "\nApproximation error: ", pi / 3 - atan_approx(x, u))

MC approx. of arctan(x^(0.5)):  1.048289 
Approximation error:  -0.001091812

xs <- seq(from = 1, to = 2, by = 0.01)
atans <- sapply(xs, atan_approx, u = u) 

plot(xs, atans, col = "red") 
lines(xs, atan(xs), col = "black") 
legend("right", legend=c("MC arctan", "R arctan"),
col=c("red", "black"), lty=c(NA, 1), pch=c(1, NA))

Question 2

Calculate the area defined by the points \(\{(x,y):|x|^{1.6}+|y|^{1.6}<1\}\).

Focus on \(S = \{(x,y): x^{1.6}+y^{1.6}<1, \, x,y \in [0,1]\}\)

Question 2

xs <- seq(-1, 1, 0.001) 
ys <- (1 - abs(xs)^(1.6))^(1 / 1.6)

plot(xs, ys, type = "l", ylim = c(-1,1), lwd = 2)
lines(xs, -ys, lwd = 2)
x_pos <- seq(0, 1, 0.001) 
y_pos <- (1 - abs(x_pos)^(1.6))^(1 / 1.6)
polygon(c(x_pos, rev(x_pos)), c(y_pos, rep(0, length(y_pos))), col = "lightgrey")
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)

Question 2

Intuitive answer

xs <- seq(0, 1, 0.001) 
ys <- (1 - abs(xs)^(1.6))^(1 / 1.6)

plot(xs, ys, type = "l", lwd = 2)
x_pos <- seq(0, 1, 0.001) 
y_pos <- (1 - abs(x_pos)^(1.6))^(1 / 1.6)
polygon(c(x_pos, rev(x_pos)), c(y_pos, rep(0, length(y_pos))), col = "lightgrey")

set.seed(123) 
N <- 500 
xs <- runif(N)
ys <- runif(N) 
inds <- xs^1.6 + ys^1.6 < 1 
points(xs[inds], ys[inds], col = "blue")
points(xs[!inds], ys[!inds], col = "red")

\[ \Pr(U \in S) = \frac{\textrm{Area of } S}{\textrm{Area of } [0,1] \times [0,1]} = \textrm{Area of } S, \quad U \sim \mathcal{U}((0,1)^2) \]

Question 2

\[ \begin{aligned} \textrm{Area of S} &= \Pr(U \in S) \\ &= \textrm{E}\left[\mathbb{1} \left\{U_1^{1.6} + U_2^{1.6} < 1 \right\} \right] \\ &\approx \frac{1}{n} \sum_{i = 1}^n \mathbb{1}\left\{U_1^{1.6} + U_2^{1.6 }< 1 \right\}, \quad U \sim \mathcal{U}((0,1)^2) \end{aligned} \] We can simulate \(U = (U_1, U_2)\) by taking two random samples \(U_i \sim \mathcal{U}(0,1)\).

Question 2

Analytical answer

\[ \begin{aligned} \textrm{Area of S} %&= \int_0^1 (1 - x^{1.6})^{1/1.6} dx \\ &= \int_0^1 \int_0^{(1-x^1.6)^{1/1.6}} 1 dy dx \\ &= \int_0^1 \int_0^1 \mathbb{1}\left\{y < (1-x^1.6)^{1/1.6}\right\} dy dx \\ &= \int_0^1 \int_0^1 \mathbb{1}\left\{x^{1.6}+y^{1.6}<1\right\} dy dx \\ &= \textrm{E}\left[\mathbb{1} \left\{U_1^{1.6} + U_2^{1.6 }< 1 \right\} \right] \\ \end{aligned} \]

N <- 10^5 
set.seed(123) 
X <- runif(N) 
Y <- runif(N) 
Z <- X^1.6 + Y^1.6 < 1 
pr <- mean(Z) 

area <- 4 * pr 
cat("Area: ", area) 

Area:  2.83872

Question 3

Use the function hist, boxplot, qqnorm, qqline to plot a random variable \(z\) from standard normal distribution and from a normal distribution with mean \(20\) and standard deviation \(5\) of sample size \(10^3\). More specifically, plot hist(z), plot(1:1000,z), qqnorm(z).

par(mfrow = c(2, 2), asp = 1)
z <- rnorm(10^3, mean = 0, sd = 1)
hist(z, col = "steelblue3")
qqnorm(z, pch = 1, cex = 1, col = "mediumpurple1")
qqline(z)
boxplot(z, main = "Box Plot", col = "tomato")
plot(1:10^3, z, xlab = "N", ylab = "z", type = "l", main = "A white noise process") 

Question 4

1. Create a random sample of size \(1000\) from a normal distribution with mean \(10\) and standard deviation \(20\). Compute the mean, the standard deviation, the median, and the quantiles at \(5\%\) and \(95\%\) probability level.
2. Write a function that returns the moving averages: \[ \frac{x_1 + x_2 + x_3}{3}, \, \frac{x_2 + x_3 + x_4}{3}, \, \ldots, \, \frac{x_{n-2} + x_{n-1} + x_{n}}{3} \,. \] Use it for the sample in (a) to get a histogram of values.

Remark: Moving Average processes are just averages of white noise.

Question 4.a

set.seed(123)
X <- rnorm(1000, mean = 10, sd = 20)
m <- mean(X) 
std <- sd(X) 
med <- median(X) 
qs <- c(0.05, 0.95) 
quants <- quantile(X, probs = qs) 
cat("Mean: ", m, "sd: ", std, "Median: ", med, "Quantiles: ", quants) 

Mean:  10.32256 sd:  19.8339 Median:  10.18419 Quantiles:  -22.45169 43.52268

Question 4.b

Vectorise the computation of the moving average:

\[ \textrm{MA}_3 = \frac{1}{3} \left( \begin{bmatrix} x_1 \\ \vdots \\ x_{n-2} \end{bmatrix} + \begin{bmatrix} x_2 \\ \vdots \\ x_{n-1} \end{bmatrix} + \begin{bmatrix} x_3 \\ \vdots \\ x_{n} \end{bmatrix} \right) \]

MA3_f <- function(xVector){
    n <- length(xVector) 
    x <- (xVector[-c(n-1, n)] + xVector[-c(1, n)] + xVector[-c(1, 2)]) / 3 
    return (x)
}
MA <- MA3_f(X)
hist(MA,col = "steelblue3", main = "Histogram of Moving Average")

Question 5

Let \(Y_t\) be the yield of a fund at \(t\) year and \(Y_t = \exp(X_t)\), where \(X_1, X_2, \ldots , X_t\) is an independent identically distributed sequence of normal random variables with mean \(0.09\) and variance \(0.1\).

An investor invests an amount \(1\) at times \(t = 0, 1, 2\) and wants to know what is the probability that the accumulation at \(t = 3\) will be higher than \(4\) for a vector of \(10^6\) possible returns.

\[ \textrm{Accumulation at } t_3 \, = \exp\{X_3\} + \exp\{X_2 + X_3\} + \exp\{X_1 + X_2 + X_3\} \]

prb <- function(threshold, X = NULL, N = 10^6, seed = 123){
    if(is.null(X)){
        set.seed(seed)
        X1 <- rnorm(N, 0.09, sqrt(0.1)) 
        X2 <- rnorm(N, 0.09, sqrt(0.1))
        X3 <- rnorm(N, 0.09, sqrt(0.1)) 
        X <- exp(X1 + X2 + X3) + exp(X2 + X3) + exp(X3)
    }
  mean(X > threshold)
}

cat("Probability that accumulation is > 4 at t = 3: ", prb(4)) 

Probability that accumulation is > 4 at t = 3:  0.411847