Financial Mathematics 1 - Excercise 1

Exercise 1.1

This exercise guides you through the basics of R.

a) First steps. Type in the following commands in the console (end each line by hitting ENTER) and see what happens.

b) Matrices, loops, writing and reading tables. Create two 3 × 3 matrices, x and y, carrying values of your choice. Type help(dim), help(matrix) and help(array) to find out how to create matrices. Check the result of the naive multiplication x * y. Write a loop that properly multiplies the two matrices. Find out whether R has built-in support for matrix multiplication (without using your loop). Use write.table to save your result matrix to the hard drive and read.table to read it back in from there.

library("tidyverse")

x <- matrix(1:9, nrow = 3, ncol = 3)
y <- matrix(10:18, nrow = 3, ncol = 3)
x
y
x * y


# Initialize a new matrix C with the same dimensions as x and y
C <- matrix(0, nrow = nrow(x), ncol = ncol(y))
# Nested loop to multiply each element of x by the corresponding element of y and store in C
for (i in 1:nrow(x)) {
  for (j in 1:ncol(y)) {
    for (k in 1:ncol(x)) {
      C[i, j] <- C[i, j] + x[i, k] * y[k, j]
    }
  }
}
C

x %*% y

A matrix: 3 × 3 of type int

1	4	7
2	5	8
3	6	9

A matrix: 3 × 3 of type int

10	13	16
11	14	17
12	15	18

A matrix: 3 × 3 of type int

10	52	112
22	70	136
36	90	162

A matrix: 3 × 3 of type dbl

138	174	210
171	216	261
204	258	312

A matrix: 3 × 3 of type dbl

138	174	210
171	216	261
204	258	312

c) # ===== Cell 4: [code] ===== f <- function(x,y){ cos(y)/(1+x^2)}

f <- function(x,y) {
    cos(y)/(1+x^2)
}
matrix_mult <- function(n,x,y){
  # Initialize a new matrix C with the same dimensions as x and y
  C <- matrix(0, nrow = nrow(x), ncol = ncol(y))
  # Nested loop to multiply each element of x by the corresponding element of y and store in C
  for (i in 1:nrow(x)) {
    for (j in 1:ncol(y)) {
      for (k in 1:ncol(x)) {
        C[i, j] <- C[i, j] + x[i, k] * y[k, j]
      }
    }
  }
  return(C)
}

d) Random numbers. Type help(rnorm) to learn how to create normally distributed random numbers. Compute the mean and the standard deviation of a vector x of normally distributed random numbers. Consider also the histogram plot, hist(x). Repeat these computations for random vectors of other distributions.

# help(rnorm)
mu <- 0
sigma <- 1
n <- 100
x <- rnorm(n, mean = mu, sd = sigma)
poiss <- rpois(n, lambda = 1)
cauchy <- rcauchy(n, location = 0, scale = 1)
expo <- rexp(n, rate = 1)

hist(x)
hist(poiss)
hist(cauchy)
hist(expo)

e) If-clauses. The if-clause allows you to dynamically react to the current value of a variable. ===== Cell 6: [code] ===== x <- runif(1,-1,1) if(x>=0){ print("x nonnegative") }else{ print("x negative") } ===== Cell 7: [markdown] ===== Write a program that loops through a vector and replaces negative numbers by their absolute value using if-clauses.

n <- 10
x <-rnorm(n, 0, 1)
y <- 1:n
for (i in 1:length(x)) {
    y[i] <- abs(x[i])
}
x
y

1. 0.98797039299508
2. -1.32012765890516
3. 1.32985233471214
4. -0.149971608388607
5. 0.113215685407754
6. 0.658643994680162
7. 0.965391908632048
8. -1.49347080632981
9. -1.32862056538383
10. 0.337193400965796

1. 0.98797039299508
2. 1.32012765890516
3. 1.32985233471214
4. 0.149971608388607
5. 0.113215685407754
6. 0.658643994680162
7. 0.965391908632048
8. 1.49347080632981
9. 1.32862056538383
10. 0.337193400965796

f) Plots. R offers several ways to plot data and results. # ===== Cell 8: [code] ===== gridpoints <- 100 xgrid <- seq(-10,10,20/gridpoints) ypoints <- cos(axis) plot(axis, ypoints)

axis <- seq(-2,2,0.01)
gridpoints <- 100
xgrid <- seq(-10,10,20/gridpoints)
ypoints <- cos(axis)
plot(axis, ypoints)

Exercise 1.2

Write functions in R that integrate a function f : [0, 1] →[0, 1] using
    1. the trapezoidal rule 
    2. Monte Carlo integration 
    3. The R-command integrate Test your functions using f(x) = | sin(x)| and compare the results.

Trapezoidal rule:

To approximate $\int_a^b f(x)dx$, replace $f$ by a piecewise linear interpolant and integrate that instead. Each piece's integral is the area of a trapezoid. Split $[a,b]$ into $n$ subintervals of equal widt: $$ h = \frac{b-a}{n}, \ \ \ x_i = a+ih, \ \ \ i = 0, \dots n $$ $$ \int_a^b f(x) dx \simeq T_n = \frac{h}{2} \left[ f(x_o) + 2 f(x_1) + \dots + 2 f(x_{n-1}) + f(x_n) \right] $$

Monte Carlo Integration

Monte Carlo integration treats an integral as an expected value. For function $f$ on $[0,1]$ with $U \sim \text{Uniform}(0,1)$: $$ \int_0^1f(x)dx=\mathbb{E}[f(U)] $$ Averaging $f$ at random points estimates that average. Draw $U_1,\dots,U_n$ i.i.d. uniform on $[0,1]$ and set $$ \hat{I}_n = \frac{1}{n}\sum_{i=1}^nf(U_i) $$ By the law of large numbers $\hat{I}_n \to \int_0^1 f(x)dx$ almost surely.

trap <- function(f, a, b, n) {
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  integral <- (h / 2) * (y[1] + 2 * sum(y[2:n]) + y[n + 1])
  return(integral)
}

monte_carlo <- function(f, a, b, n) {
  x <- runif(n, min = a, max = b)
  integral <- (b - a) * mean(f(x))
  return(integral)
}

set.seed(1)
f <- function(x) {abs(sin(x))}
n = 100
trap(f,0,1,n)
monte_carlo(f,0,1,n)
integrate(f,0,1)

0.459693863311358

0.47788344660118

0.4596977 with absolute error < 5.1e-15

Excercise 1.3

Consider two stocks whose annual returns R1 and R2 are independent and N(µ, σ2)-distributed with µ = 0.06 and σ = 0.25, i.e., each stock has an expected return of 6% per year and the return has a standard deviation of 25%. Call set.seed(1) at the beginning of your script.

a) Simulate n = 10 000 realizations of (R1, R2) and store them in an n × 2 matrix R (one column per stock). Compute the mean and the standard deviation of each column. Hint: apply(R, 2, sd) applies sd to every column of R.

##  a)
set.seed(1)
# N(0.06, 0.25^2)
mu = 0.06
sig = 0.25
n = 10000

R <- matrix(data = rnorm(2*n, mu, sig), ncol = 2)
mu_z = apply(R, 2, mean)
sig_z = apply(R, 2, sd)
mu_z
sig_z

1. 0.058365740134584
2. 0.0589524835782514

1. 0.253089113176178
2. 0.247694419809642

b) An investor puts the fraction w of her wealth into stock 1 and 1−w into stock 2, for different weights w ∈{0, 0.1, . . . , 1}. Given a weight w, the portfolio return (consisting only of the two stocks) is R(w) = wR1 +(1−w)R2. Store all 11 possible weight vectors (w, 1−w)⊤as the columns of a 2×11 matrix W and compute the returns of all 11 portfolios (this is possible with a single matrix product %*%). Compute the mean and the standard deviation of each portfolio return and plot the standard deviation against w. Which portfolio has the lowest risk? Does splitting the investment reduce the expected return?

w <- matrix(c(seq(0,1,0.1), 1- seq(0,1,0.1)), ncol = 2)
w <- t(w)
R_w = R %*% w
colnames = rep("", 11)
for (i in 1:ncol(w)) {
    colnames[i] = paste(w[1,i], w[2,i], sep = ",")
}
colnames(R_w) = colnames
head(w)
head(R_w)

par = matrix(c(apply(R_w,2, mean), apply(X = R_w, 2, sd)), ncol = 2)
colnames(par) = c("mu", "sigma")

plot(w[1,], par[,1], type = "b", main = "Expected return vs Portfolio weight")
plot(w[1,], par[,2], type = "b", main = "Volatility vs Portfolio weight")

A matrix: 2 × 11 of type dbl

0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
1	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0

A matrix: 6 × 11 of type dbl

0,1	0.1,0.9	0.2,0.8	0.3,0.7	0.4,0.6	0.5,0.5	0.6,0.4	0.7,0.3	0.8,0.2	0.9,0.1	1,0
-0.14108290	-0.13663596	-0.13218901	-0.127742066	-0.12329512	-0.11884818	-0.11440123	-0.10995429	-0.10550734	-0.10106040	-0.09661345
-0.20413141	-0.17312719	-0.14212296	-0.111118740	-0.08011452	-0.04911029	-0.01810607	0.01289816	0.04390238	0.07490661	0.10591083
-0.19884895	-0.19385477	-0.18886059	-0.183866408	-0.17887223	-0.17387805	-0.16888387	-0.16388969	-0.15889551	-0.15390133	-0.14890715
-0.23639009	-0.16686906	-0.09734803	-0.027827002	0.04169403	0.11121506	0.18073608	0.25025711	0.31977814	0.38929917	0.45882020
-0.06510988	-0.04436120	-0.02361251	-0.002863832	0.01788485	0.03863353	0.05938221	0.08013090	0.10087958	0.12162826	0.14237694
-0.07124718	-0.07863417	-0.08602116	-0.093408153	-0.10079514	-0.10818214	-0.11556913	-0.12295612	-0.13034311	-0.13773010	-0.14511710

c) Write your own function dens(x, m, s) for the density of the N(m, s2) distribution (without using dnorm). Assuming normally distributed returns with the sample mean and sample standard deviation from b), use integrate to compute the probability of losing more than 30% in one year, P(R(w) < −0.3), for stock 1 alone (w = 1) and for the 50/50 portfolio (w = 0.5). What do you observe? Compare your results with pnorm and with the empirical frequency of such losses in your simulated data.

dens <- function(x, m, s) {
    return(1/(s*sqrt(2*pi)) * exp(-0.5*((x-m)/s)^2))
}

# P(R < -0,3)
integrate(dens, lower = -Inf, upper = -0.3, m = mu, s = sig)
pnorm(-0.3, mean = mu, sd = sig, lower.tail = TRUE)

sum(R_w[, "0.5,0.5"] < -0.3)/ nrow(R_w)
pnorm(-0.3, mean = mu, sd = sqrt((sig^2)/2), lower.tail = TRUE)

f_dens <- function(x, m, s, delta) {
    if (x > delta) {
        nominator = dens((x-m)/s, m, s)
        arg = (delta - m)/s
        denominator = s * (1 - integrate(dens, lower = -Inf, upper = arg, m=m, s=s)$value)
        return (nominator/denominator)
    }
    return (0)
}

0.0749337 with absolute error < 6.9e-08

0.0749336995343271

0.0213

0.0208517151997176

d) Write a function, rcnorm, with arguments n, mean, sd, and delta, that draws n iid samples from the distribution with the density f(x) = σ · (1 −Φ([δ −µ]/σ))ϕ([x −µ]/σ)1{x>δ}, where ϕ and Φ are the density and cumulative distribution function of a standard normal distri- bution, and µ, δ ∈R, σ > 0 are free parameters. Try to find an intuition for the distribution with density f. Can you find an example where this distribution might appear within the context of this exercise? Hint: To create one sample, use an infinite while-loop in which a normal sample is drawn; break out of the while-loop if the drawn number is larger than δ. © 2026 Technical University of Munich, Chair of Mathematical Finance.

rcnorm <- function(n, m, s, delta) {
    grid = seq(m-10*s, m+10*s, length.out = 10001)
    Fgrid <- cumsum(sapply(grid, f_dens, m = m, s = s, delta = delta))
    Fgrid <- Fgrid / Fgrid[10001]

    u <- runif(n = n, min = 0, max = 1)
    rc <- approx(x = Fgrid, y = grid, xout = u)$y
    return (rc)
}

y <- rcnorm(n = 10000, m = mu, s = sig, delta = -0.3)
y2 <- rnorm(n = 10000, mean = mu, sd = sig)
hist(y)
hist(y2)

Warning message in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
"Reduktion auf einmalige 'x' Werte"