# ST303 Stochastic Simumlation
Class 2

Philipp Sterzinger

p.sterzinger@lse.ac.uk
Office hours: Thursday, COL 8.05

16 October 2025

## 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

Code

    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

Code

    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

Code

    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

Code

    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} \]

Code

    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)`.

Code

    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

Code

    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) \]

Code

    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\} \]

Code

    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