Finanční matematika pomocí R
Jan Vávra
2025-09-23
Quandl
Máte-li účet na Nasdaq můžete
si pomocí knihovny Quandl stahovat data cen nejrůznějších
komodit
library("Quandl")
Quandl.api_key("dpruBx3AnT7Nys1ezz2N")
mydata <- Quandl("LBMA/GOLD", start_date="2015-12-31", end_date="2022-06-30")
names(mydata) = c("datum","USD_AM","USD_PM","GBP_AM","GBP_PM","EURO_AM","EURO_PM")
plot(datum, USD_AM, data = mydata, type="l",
xlab = "Rok", ylab="USD_AM", main="Vývoj otevíracích cen zlata v USD")
library(dplyr)
mydata <- mutate(mydata,
return = (USD_AM - lag(USD_AM, order_by = datum))/lag(USD_AM, order_by = datum))
head(mydata)
### Úloha
# 1) Odstraňte všechny sloupce až na "datum", "USD_AM" and "return".
# 2) Odstraňte řádky s příliš odlehlými hodnotami return, tedy |return| > 0.05.
# 3) Vykreslete histogram returnů.
par(mfrow=c(1,2), mar = c(4,4,2.5,0.5))
mydata <- mydata[, -c(3:7)]
hist(mydata$return)
# tohle nás zbaví extrémních hodnot
mydataCut <- mydata[-which(mydata$return< -0.05 | mydata$return > 0.05), ]
# nebo alternativně
mydataCut <- mydata[abs(mydata$return) <= 0.05, ]
hist(mydataCut$return)VaR, Expected Shortfall
Na stránce kaggle.com si můžete stáhnout historická data různých indexů, třeba S&P500.
remove(list=ls())
data_dir <- file.path(dirname(getwd()), "data")
# Stažení dat z mých stránek:
download.file("https://www.karlin.mff.cuni.cz/~vavraj/Rko/data/SP500.csv",
file.path(data_dir, "SP500.csv"))
mydata <- read.csv(file.path(data_dir, "SP500.csv"))
dim(mydata)## [1] 23323 7head(mydata)## Date Open High Low Close Adj.Close Volume
## 1 1927-12-30 17.66 17.66 17.66 17.66 17.66 0
## 2 1928-01-03 17.76 17.76 17.76 17.76 17.76 0
## 3 1928-01-04 17.72 17.72 17.72 17.72 17.72 0
## 4 1928-01-05 17.55 17.55 17.55 17.55 17.55 0
## 5 1928-01-06 17.66 17.66 17.66 17.66 17.66 0
## 6 1928-01-09 17.50 17.50 17.50 17.50 17.50 0summary(mydata)## Date Open High Low Close Adj.Close
## Length:23323 Min. : 4.40 Min. : 4.40 Min. : 4.40 Min. : 4.40 Min. : 4.40
## Class :character 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86
## Mode :character Median : 99.50 Median : 100.31 Median : 98.72 Median : 99.50 Median : 99.50
## Mean : 486.82 Mean : 489.69 Mean : 483.76 Mean : 486.92 Mean : 486.92
## 3rd Qu.: 834.03 3rd Qu.: 843.75 3rd Qu.: 822.89 3rd Qu.: 834.71 3rd Qu.: 834.71
## Max. :3564.74 Max. :3588.11 Max. :3535.23 Max. :3580.84 Max. :3580.84
## Volume
## Min. :0.000e+00
## 1st Qu.:1.280e+06
## Median :1.734e+07
## Mean :7.692e+08
## 3rd Qu.:5.379e+08
## Max. :1.146e+10# Výpočet denních returnů
mydata$return <- c(NA, diff(mydata$Open)) / c(NA, mydata$Open[-length(mydata$Open)])
mydata <- mydata[complete.cases(mydata), ]Historicky se Value at Risk na úrovni pravděpodobnosti \(p\) počítá jako
- \(p\)-kvantil záporných returnů či jako
- záporná hodnota \((1-p)\)-kvantil returnů.
Moderně se k výpočtu používá předpokladu normálního rozdělení,
případně zobecnění na nenormální případ, vizte
help(VaR):
library(PerformanceAnalytics)
VaR(mydata$return, p=0.95, method="historical") # historická metoda použití dat## [,1]
## VaR -0.01666246quantile(mydata$return, probs=c(0.05,0.5,0.95))## 5% 50% 95%
## -0.0166624622 0.0005290903 0.0163778699VaR(mydata$return, p=0.95, method="gaussian") # metoda používající kvantily normálního rozdělení## [,1]
## VaR -0.01907882qnorm(0.05, mean(mydata$return), sd(mydata$return))## [1] -0.01907923VaR(mydata$return, p=0.95, method="modified") # metoda používající kvantil odhadnutého rozdělení## [,1]
## VaR -0.01529831print(c(VaR(mydata$return, p=0.95, method="historical"),
VaR(mydata$return, p=0.95, method="gaussian"),
VaR(mydata$return, p=0.95, method="modified")))## [1] -0.01666246 -0.01907882 -0.01529831Expected Shortfall na pravděpodobnosti \(p\) je záporná hodnota očekávané hodnoty returnu, pakliže je return menší než jeho \((1-p)\)-kvantil.
# Průměr 5 % nejnižších pozorovaných hodnot - historická metoda
mean(mydata$return[which(mydata$return < quantile(mydata$return, 0.05))])## [1] -0.02815835ES(mydata$return, method="historical") # funkce počítající to samé## [,1]
## ES -0.02815835ES(mydata$return, method="gaussian") # funkce vycházející z normálního rozdělení## [,1]
## ES -0.02400059ES(mydata$return, method="modified") # pro nenormální rozdělení - defaultní hodnota## [,1]
## ES -0.01529831Drobná úloha na procvičení
- Nasimulujte 555 pozorování ztrát z
- exponenciálního rozdělení s parametrem
rate = 1, - exponenciálního rozdělení s parametrem
rate = 20.
- Vykreslete histogramy těchto dvou výběrů vedle sebe.
- Spočtěte historické VaR a ES pro tyto dva výběry a zamyslete se proč
Rkové funkce
VaRaESzde nefungují.
set.seed(123456789)
sample_exp_1 = rexp(555, rate=1)
sample_exp_2 = rexp(555, rate=20)
par(mfrow=c(1,2), mar = c(4,4,2.5,0.5))
hist(sample_exp_1, main="Exp(1)")
hist(sample_exp_2, main="Exp(20)")
VaR(sample_exp_1, p=0.95, method="historical") # nefunguje, kladné hodnoty VaR## VaR calculation produces unreliable result (inverse risk) for column: 1 : -0.047686057648395## [,1]
## VaR NAVaR(-sample_exp_1, p=0.95, method="historical") # nefunguje, kvantil pod -1 ## VaR calculation produces unreliable result (risk over 100%) for column: 1 : 3.01853769032245## [,1]
## VaR -1-quantile(sample_exp_1, 0.95) # jedná se o vektor ZTRÁT --> 0.95## 95%
## -3.018538-quantile(sample_exp_2, 0.95)## 95%
## -0.1457132mean(sample_exp_1[which(sample_exp_1 > quantile(sample_exp_1, 0.95))])## [1] 4.31829mean(sample_exp_2[which(sample_exp_2 > quantile(sample_exp_2, 0.95))])## [1] 0.1979028Internal Rate of Return
Takto pro danou posloupnost cash flows v Rku spočteme internal rate of return:
remove(list=ls())
library(FinancialMath)
# ?IRR
# IRR # tělo funkce
IRR(cf0=100, # cash flow v nultém čase
cf=c(10,10,110), # ostatní cash flow
times=1:3, # časy, kdy k nim má dojít
plot = FALSE) # zda má vykreslit diagram## [1] 0.1IRR(cf0=100, cf=c(20,10,110), times=1:3, plot = FALSE)## [1] 0.1377142Používá funkci polyroot pro výpočet kořenu polynomiální
rovnice.
Net Present Value při dané úrokové míře \(i\):
# ?NPV
par(mfrow=c(1,1))
NPV(cf0=100, cf=c(10,10,110), times=1:3, i=0.11, TRUE)
## [1] -2.443715Ukažme si na delší časové řadě cash-flows, jak by se to dělalo ručně:
cashflows <- data.frame(time = 0:5,
cf = c(-100,5,20,25,30,50),
interest_rate = rep(0.05,6))
print(cashflows)## time cf interest_rate
## 1 0 -100 0.05
## 2 1 5 0.05
## 3 2 20 0.05
## 4 3 25 0.05
## 5 4 30 0.05
## 6 5 50 0.05# ruční výpočet NPV
cashflows <- transform(cashflows, pv = cf / (1+interest_rate)^time)
npv <- sum(cashflows$pv)
round(npv,5)## [1] 8.35582NPV(cf0=cashflows$cf[1],
cf=cashflows$cf[-1],
times=cashflows$time[-1],
i=0.05)## [1] 8.355817# IRR = Taková úroková míra, při níž je NPV nulové.
# Zde už by to byl polynom 5.řádu, což ručně neupočítáme.
# Lze použít polyroot, jen je těžší najít koeficienty polynomu.
# Ruční výpočet pomocí uniroot:
irr <- uniroot(function(irr){sum(cashflows$cf / (1+irr)^{cashflows$time})},
interval = c(0,1))
IRR(cf0=cashflows$cf[1],
cf=cashflows$cf[-1],
times=cashflows$time[-1])## [1] 0.07315677# Znázornění na obrázku
rates <- seq(0, 0.5, 0.01) # hustá síť různých úrokových měr
npvs <- sapply(rates, function(i){sum(cashflows$cf / (1+i)^{cashflows$time})})
par(mfrow = c(1,1), mar = c(4,4,2.5,0.5))
plot(npvs~rates, type="l")
abline(h=0, col="blue")
abline(v=c(0,0.05,0.1,0.15,0.2), lty=3)
abline(v=irr$root, lty=5, lwd=2)
Další úloha uvažuje různé úrokové míry v různých časech, jaká je NPV?
Hint: Použijte funkci cumprod.
cashflows$interest_rate <- c(0.0, # aby (1/(1+i)) bylo 1
0.25,0.25,0.15,0.15,0.125) # úrokové míry v dalších obdobích
cashflows$discount_factor <- cumprod(1/(1+cashflows$interest_rate))
cashflows## time cf interest_rate pv discount_factor
## 1 0 -100 0.000 -100.000000 1.0000000
## 2 1 5 0.250 4.761905 0.8000000
## 3 2 20 0.250 18.140590 0.6400000
## 4 3 25 0.150 21.595940 0.5565217
## 5 4 30 0.150 24.681074 0.4839319
## 6 5 50 0.125 39.176308 0.4301617cashflows <- transform(cashflows, pv = cf * discount_factor) # diskontované cash-flows
(npv2 <- sum(cashflows$pv)) ## [1] -33.26091Náhodné procesy a simulace
Opět si vezměme index S&P500.
remove(list=ls())
data_dir <- file.path(dirname(getwd()), "data")
mydata <- read.csv(file.path(data_dir, "SP500.csv"))
head(mydata)## Date Open High Low Close Adj.Close Volume
## 1 1927-12-30 17.66 17.66 17.66 17.66 17.66 0
## 2 1928-01-03 17.76 17.76 17.76 17.76 17.76 0
## 3 1928-01-04 17.72 17.72 17.72 17.72 17.72 0
## 4 1928-01-05 17.55 17.55 17.55 17.55 17.55 0
## 5 1928-01-06 17.66 17.66 17.66 17.66 17.66 0
## 6 1928-01-09 17.50 17.50 17.50 17.50 17.50 0summary(mydata)## Date Open High Low Close Adj.Close
## Length:23323 Min. : 4.40 Min. : 4.40 Min. : 4.40 Min. : 4.40 Min. : 4.40
## Class :character 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86 1st Qu.: 23.86
## Mode :character Median : 99.50 Median : 100.31 Median : 98.72 Median : 99.50 Median : 99.50
## Mean : 486.82 Mean : 489.69 Mean : 483.76 Mean : 486.92 Mean : 486.92
## 3rd Qu.: 834.03 3rd Qu.: 843.75 3rd Qu.: 822.89 3rd Qu.: 834.71 3rd Qu.: 834.71
## Max. :3564.74 Max. :3588.11 Max. :3535.23 Max. :3580.84 Max. :3580.84
## Volume
## Min. :0.000e+00
## 1st Qu.:1.280e+06
## Median :1.734e+07
## Mean :7.692e+08
## 3rd Qu.:5.379e+08
## Max. :1.146e+10Zkusíme najít vhodný model pro tento index. Abychom byli schopni vyhodnotit jeho kvalitu, tak si nechme posledních 150 hodnot stranou:
training <- head(mydata, length(mydata[,1])-150)
testing <- tail(mydata, 150)Budeme jej modelovat jako Geometric Brownina Motion, tedy předpokládáme, že denní log-returny jsou iid (nezávislé stejně rozdělené) z normálního rozdělení. Nejprve je tedy zapotřebí spočíst returny a log-returny \(\log(1+\mathsf{return})\):
# return_t = (price_t - price_{t-1}) / price_{t-1}
# return_t = price_t / price_{t-1} -1
returns <- diff(training$Open) / training$Open[-length(training$Open)]
# logret_t = log(1+return_t) = log(price_t / price_{t-1})
logret <- log(1+returns)Nyní odhadněme průměr a směrodatnou odchylku z trénovacího datasetu:
mu <- mean(logret)
sigma <- sd(logret)Je náš model rozumnou volbou pro tato data? Pomožme si QQ-plotem, kde napozorované kvantily jsou vykresleny proti těm teoretickým za předpokladu normálního rozdělení. Za splnění předpokladu by tedy měly ležet ideálně v jedné přímce:
par(mfrow = c(1,1), mar = c(4,4,2.5,0.5))
qqnorm(logret)
qqline(logret, col = 2)
Generujme si nových 150 hodnot z predikovaného rozdělení:
set.seed(123456789)
logReturns <- rnorm(150, mu, sigma)A nyní se vraťme k původní řadě a nagenerujme vývoj na dalších 150 dní:
actPrice <- training[dim(training)[1], "Open"] # poslední známá cena v trénovacím datasetu
# exp(logret_t) = price_t / price_{t-1} --> price_t = price_{t-1} * exp(logret_t)
simulatedSP <- cumprod(c(actPrice, exp(logReturns)))Porovnejme, jak se naše náhodně vygenerovaná řada liší od skutečnosti:
par(mfrow=c(2,1), mar=c(4,4,0.5,0.5))
plot(testing$Open, type="l", col="red", ylim = range(c(testing$Open, simulatedSP)))
lines(simulatedSP, type="l")
fullTrajectory <- c(training$Open, testing$Open)
fullTrajectory_simul <- c(training$Open, simulatedSP)
plot(fullTrajectory, type="l", col="red", ylim = range(fullTrajectory, fullTrajectory_simul))
lines(fullTrajectory_simul, type="l")
Nagenerujte 100 různých variant predikcí a porovnejte s testovacím datasetem. Napočítejte 95% predikční interval pro poslední známou hodnotu.
ruzne_scenare <- replicate(100, cumprod(c(actPrice, exp(rnorm(150,mu,sigma)))))
ruzne_scenare <- ruzne_scenare[-1,]
# Napočítání dolních a horních kvantilů - pro každý časový okamžik
predikcni_intervaly <- apply(ruzne_scenare, 1, quantile, probs = c(0.025, 0.975))
# Vykreslení jednotlivých scénářů
par(mfrow=c(1,1), mar=c(4,4,0.5,0.5))
plot(testing$Open, type="l", col="red", ylim = range(ruzne_scenare), lwd = 2) # reálný scénář
for(i in 1:100){
lines(ruzne_scenare[,i], type="l", lwd = 0.5) # jednotlivé generované scénáře
}
# Zakreslení predikčních intervalů jako polygon s průsvitnou barvou (alpha = 0.5)
polygon(x = c(1:150, 150:1),
y = c(predikcni_intervaly[1,], rev(predikcni_intervaly[2,])),
border = "blue", col = rgb(0,0,1,alpha=0.5))
Více o analýze a simulování stock market returnů si můžete přečíst zde
či zde.
Dodatečné úlohy k úročení
Úloha 1
Vítěz loterie má často možnost vybrat si, zda chce peníze dostávat v rovnoměrných splátkách nebo jako jednorázovou částku. Předpokládejme, že vyhrajete v loterii 50 000 dolarů a máte dvě možnosti:
- Pět rovnoměrných splátek po 10 000 dolarech v časech 0,…,4;
- Jednorázová částka 45 000 dolarů okamžitě.
Předpokládejme, že banka vyplácí 5% úrokovou sazbu z vkladu za dané období. Která z možností 1) a 2) je výhodnější? Změnila by se odpověď, pokud by pět splátek začalo v čase 1?
pmt = 10000
n = 5
i = 0.05
cf1 <- data.frame(time=(0:4),
pmt = pmt)
cf1$pv <- cf1$pmt/(1+i)^cf1$time
sum(cf1$pv)## [1] 45459.51cf2 <- data.frame(time=(1:5),
pmt = pmt)
cf2$pv <- cf2$pmt/(1+i)^cf2$time
sum(cf2$pv)## [1] 43294.77# or using explicit formula for PV of annuity (or annuity due)
pv1 <- pmt*(1+i)*(1-(1+i)^(-5))/i
pv2 <- pmt*(1-(1+i)^(-5))/i
pv1## [1] 45459.51pv2## [1] 43294.77# or using NPV:
library(FinancialMath)
npv1 <- -NPV(10000, rep(-10000,4), times=1:4, i=0.05)
npv2 <- NPV(0, rep(10000,5), times=1:5, i=0.05)Úloha 2
Americká ropná a plynárenská společnost „Pioneer Natural Resources Co. (PXD)“, která patří mezi společnosti indexu S&P500, má vyplácet roční dividendu ve výši 3,89 USD na akcii a společnost bude tyto dividendy vyplácet na dobu neurčitou.
Kolik jsou investoři ochotni zaplatit za akcii s požadovanou mírou výnosu 6 %?
- dividenda je vyplácena ročně,
- dividenda je vyplácena čtvrtletně – bylo by nutné přepočítat nominální úrokovou sazbu, (předpokládali bychom, že požadovaný výnos 6 % je efektivní úroková sazba).
pmt <- 3.89
i <- 0.06
(price <- pmt/i)## [1] 64.83333pmtQuart <- pmt/4
iQuart <- (1+i)^(1/4)-1
(priceQuart <- pmtQuart/iQuart) #NPV for perpetuity## [1] 66.27438Úloha 3
Předpokládejme, že vlastníte bytový komplex. Nájemníci vám každý rok platí pevné nájemné. Nájemné můžete zvýšit, abyste vyrovnali vyšší inflaci. Tento peněžní tok lze tedy popsat jako anuitní platbu, která roste tempem inflace.
Vypočítejte současnou a budoucí hodnotu této anuity po 5 letech za předpokladu, že
- nájemníci platí nájemné pozadu (polhůtné),
- celková částka, kterou nájemníci zaplatí v tomto roce, je 1000 $,
- očekáváte, že míra inflace bude každoročně na úrovni 3 %,
- úroková sazba je 0,06 (= 6 %).
pmt <- 1000
growth <- 0.03
n <- 5
i <- 0.06
cfs <- data.frame(time = 1:5,
cf = pmt)
cfs <- transform(cfs, cf = pmt*(1+growth)^(time))
(cfs <- transform(cfs, pv = cf/(1+i)^time))## time cf pv
## 1 1 1030.000 971.6981
## 2 2 1060.900 944.1972
## 3 3 1092.727 917.4747
## 4 4 1125.509 891.5084
## 5 5 1159.274 866.2770(npv <- sum(cfs$pv))## [1] 4591.155cfs <- transform(cfs, fv = cf*(1+i)^(5-time))
(fv <- sum(cfs$fv))## [1] 6144.002# or equivalently
i_net <- (1+i)/(1+growth) -1
npv_2 <- pmt*(1-(1+i_net)^(-5))/i_net
fv_2 <- npv_2*(1+i)^5Úloha 4
Pan Hornet zemřel a všechny jeho peníze, 2 500 000 dolarů, byly převedeny do fondu s roční nominální úrokovou sazbou i = 7 % a měsíčním úročením. Přál si, aby každé z jeho tří dětí dostalo stejnou nominální částku peněz k 18. narozeninám, takže poté, co nejmladší dítě dostane své peníze, bude fond prázdný.
Kolik peněz každé dítě obdrží k 18. narozeninám, pokud je jim v okamžiku úmrtí pana Horneta 10, 15 a 16 let? Předpokládejme, že všechny narozeniny se shodují s dnem úmrtí.
i = 0.07
m = 12
im = i/m
years_to_18 = 18 - c(10,15,16)
pv <- function(x){
-2500000 + x * sum(1/((1+im)^(m*years_to_18))) # x je stejná částka k vyplacení
}
# Kořen nalezneme buď pomocí funkce uniroot
(root <- uniroot(pv, c(700000, 1500000)))## $root
## [1] 1109666
##
## $f.root
## [1] 4.656613e-10
##
## $iter
## [1] 2
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103562e-05# nebo si uvědomíme, že jde o lineární funkci
2500000 / sum(1 / ((1+im)^(m*years_to_18)))## [1] 1109666par(mfrow = c(1,1), mar = c(4,4,0.5,0.5))
x <- seq(700000, 1500000, 100)
plot(x, sapply(x, pv), type="l")
abline(h = 0, col="blue")
points(x = root$root, y = 0, col="red", pch = 16, cex = 1.5)
abline(v = root$root, lty = 2)
Úloha 5 - IRR a NPV profile
Předpokládejme tři peněžní toky (cash-flows) A, B, C následovně:
- A = -1000, 1000, 750
- B = -700, 900, 100
- C = 1000, -700, -700
Platby se uskutečňují v časech 0, 1, 2. Úkolem je
- sestavit profil NPV pro každý CF a zobrazit je v jediném grafu,
- najít a vynést IRR (i*) pro každý CF,
- určit, pro které \(i\) je A výnosnější než C a pro které \(i\) jsou stejné.
cfA = c(-1000, 1000, 750)
cfB = c(-700, 900, 100)
cfC = c(1000, -700, -700)
times = 0:2
pv = function(i,cf){ # present value - function
sum(cf/((1+i)^times))
}
x <- seq(0, 0.8, by = 0.001) # hustá množina úrokových měr, pro které se má napočítat
pvA <- sapply(x, pv, cf = cfA) # aplikujeme funkce pv na hodnoty x
pvB <- sapply(x, pv, cf = cfB) # přičemž specifikujeme druhý parametr
pvC <- sapply(x, pv, cf = cfC)
# Nebo použitím funkce NPV z balíčku FinancialMath
pvA2 <- sapply(x, NPV, cf0 = cfA[1], cf = cfA[-1], times=times[-1]) # je třeba udat ostatní parametry
# zbylý parametr bude proměnná
all.equal(pvA, pvA2)## [1] TRUE# Maticově to lze zapsat na dva řádky a bude to fungovat pro libovolné množství cash-flows
cfall <- matrix(c(cfA, cfB, cfC), ncol = 3, byrow = TRUE)
pvall <- apply(cfall, 1, function(cashflow){sapply(x, pv, cf = cashflow)})
COLs <- c("aquamarine", "brown", "chartreuse")
names(COLs) <- LETTERS[1:3]
par(mfrow = c(1,1), mar = c(4,4,0.5,0.5))
plot(x, pvA, type="l", ylab = "PV", xlab = "interest rate", col = COLs["A"], lwd = 2)
abline(h = 0, col = "blue")
lines(x, pvB, col = COLs["B"], lwd = 2)
lines(x, pvC, col = COLs["C"], lwd = 2)
legend("topright", LETTERS[1:3], col = COLs, lwd = 2, lty = 1)
# IRR nalezneme jako kořeny funkcí
rootsA <- uniroot(pv, interval = c(0,0.8), cf = cfA) # cf je třeba dodat jako dodatečný argument
rootsB <- uniroot(pv, interval = c(0,0.8), cf = cfB)
rootsC <- uniroot(pv, interval = c(0,0.8), cf = cfC)
# nebo pomocí funkce IRR
# zde je maticové řešení na řádek a potenciálně libovolné množství cash-flows
rootsall <- apply(cfall, 1, function(cashflow){uniroot(pv, interval = c(0,0.8), cf = cashflow)$root})
roots <- c(rootsA$root, rootsB$root, rootsC$root)
points(x = roots, y = rep(0, length(roots)),
col="red", pch = 16, cex = 1)
# Nalezením průsečíku A a C získáme úrokovou míru takovou, že
# je-li skutečná nižší, tak je A více profitabilní,
# je-li skutečná vyšší, tak je C více profitabilní.
rozdil_AC <- function(i){
pv(i, cfA) - pv(i, cfC)
}
prusecik_AC <- uniroot(rozdil_AC, interval = c(0,0.8))$root
points(x = prusecik_AC, y = pv(prusecik_AC, cfA),
col = "magenta1", pch=16)
Úloha 6
Pozorujeme celkem 50 plateb konstantní hodnoty 1000 v časech 1 až 50. Pro každé časové období existuje úroková sazba, která je náhodná a pochází z rovnoměrného rozdělení definovaného na intervalu [0,01, 0,2].
- Vygenerujte takovou posloupnost úrokových sazeb pro časy 1:50 a uložte ji do smysluplného datového rámce.
- Vypočítejte aritmetický průměr úrokových sazeb.
- Vypočítejte geometrický průměr úrokových sazeb.
- Vypočítejte PV plateb pomocí
- „skutečných“ pozorovaných úrokových sazeb, jedinečných pro každý čas,
- aritmetického průměru úrokových sazeb jako konstantní úrokové sazby,
- geometrického průměru úrokových sazeb jako konstantní úrokové sazby.
- Porovnejte výsledky z předchozího úkolu. Jsou aritmetický a/nebo geometrický průměr vhodnými charakteristikami skutečných úrokových sazeb?
- Zopakujte generování úrokových sazeb a všechny výpočty 100× a vhodně zakreslete výsledky. Jak byste spočetli skutečnou střední hodnotu NPV?
rm(list=ls())
set.seed(123456789)
n <- 50
cfs <- data.frame(time = 1:n,
pmt = 1000,
i_real = runif(n, min=0.01, max=0.2))
# 1)
(cfs$i_arith <- mean(cfs$i_real))## [1] 0.1030059# 2)
geom_by_def <- prod(cfs$i_real)^{1/length(cfs$i_real)} # skoro 0^0... numericky nestabilní
# log geom_mean = 1/n sum log(interests)
(cfs$i_geomMean <- exp(mean(log(cfs$i_real))))## [1] 0.07935707all.equal(cfs$i_geomMean[1], geom_by_def) # to samé## [1] TRUE# 3)
# a)
cfs$discount_factor <- cumprod(1/(1+cfs$i_real))
cfs <- transform(cfs, pv_real = pmt * discount_factor)
npv_real <- sum(cfs$pv_real)
# b)
cfs <- transform(cfs, pv_arith = pmt/(1+i_arith)^time)
npv_arith <- sum(cfs$pv_arith)
# c)
cfs <- transform(cfs, pv_gMean = pmt/(1+i_geomMean)^time)
npv_gMean <- sum(cfs$pv_gMean)
# d) Stejný náhodně vygenerovaný úrok po všech 50 období - přesný výpočet střední hodnoty
# Uvažme i ~ Unif(a, b)
# potom E 1/(1+i)^t = 1/(b-a) * ((1+b)^(1-t) - (1+a)^(1-t))/(1-t) pro t > 0, t!=1
# = 1/(b-a) * log((1+b)/(1+a)) pro t = 1
dis_fac_mean <- function(t, a=0.01, b=0.2){
if(a >= b){
stop("a >= b, but a<b is supposed.")
}
if(t < 0){
stop("t is negative.")
}
if(t == 0){
return(1)
}
if(t == 1){
return(1/(b-a) * log((1+b)/(1+a)))
}
return(1/(b-a) *((1+b)^(1-t) - (1+a)^(1-t))/(1-t))
}
# dis_fac_mean(-1)
sapply(c(0,0.5,1,1.5), dis_fac_mean) # vyzkoušení funkce## [1] 1.0000000 0.9521848 0.9072170 0.8649080cfs <- transform(cfs, pv_inaccurate = pmt * sapply(time, dis_fac_mean))
npv_inaccurate <- sum(cfs$pv_inaccurate)
# e) Různé náhodně vygenerované úroky - přesný výpočet střední hodnoty
dis_fac_mean2 <- function(t, a=0.01, b=0.2){
# E 1 / (1+i_1) / (1+i_2) / ... / (1+i_t) = (E 1/(1+i))^t
if(a >= b){
stop("a >= b, but a<b is supposed.")
}
if(t < 0){
stop("t is negative.")
}
if(t == 0){
return(1)
}
return((1/(b-a) * log((1+b)/(1+a)))^t)
}
sapply(c(0,0.5,1,1.5), dis_fac_mean) # vyzkoušení funkce## [1] 1.0000000 0.9521848 0.9072170 0.8649080cfs <- transform(cfs, pv_accurate = pmt * sapply(time, dis_fac_mean2))
npv_accurate <- sum(cfs$pv_accurate)
# 4)
c(npv_arith, npv_gMean, npv_real, npv_inaccurate, npv_accurate) # geomMean < arithMean ...... zaručené## [1] 9636.029 12324.477 8431.967 12772.822 9702.706# 5)
simulace <- t(replicate(100, {
n <- 50
time <- 1:n
pmt <- 1000
i_real <- runif(n, min=0.01, max=0.2)
i_arith <- mean(i_real)
i_geomMean <- exp(mean(log(i_real)))
c(sum(pmt / (1+i_arith)^time),
sum(pmt / (1+i_geomMean)^time),
sum(pmt * cumprod(1/(1+i_real)))
)
}))
dim(simulace)## [1] 100 3summary(simulace)## V1 V2 V3
## Min. : 7913 Min. : 9077 Min. : 7613
## 1st Qu.: 8919 1st Qu.:10509 1st Qu.: 8813
## Median : 9422 Median :11321 Median : 9592
## Mean : 9429 Mean :11375 Mean : 9675
## 3rd Qu.: 9850 3rd Qu.:12031 3rd Qu.:10418
## Max. :11607 Max. :14672 Max. :11883par(mfrow = c(1,1), mar = c(2.5,2.5,0.5,0.5))
boxplot(simulace, names = c("Aritmetický průměr", "Geometrický průměr", "Skutečné úroky"), type = "n")
for(i in 1:dim(simulace)[1]){
lines(1:3, simulace[i,], col = "grey", lwd = 0.2)
}
boxplot(simulace, add = T, names = NA, col = c("skyblue", "darksalmon", "gold"))
abline(h = c(npv_inaccurate, npv_accurate), col = c("grey", "red"), lwd = c(0.5,2), lty = 5)
legend("topleft", c("Jedna náh. úr. míra", "Každé období jiná náh. úr. míra"),
title = "Přesný výpočet střední hodnoty",
col = c("grey", "red"), lwd = c(0.5,2), lty = 5, cex = 0.7)