rm(list = ls()) rm(umrtnostni_tabulky) rm(j,i) ######################################################################## ## import umrtnostnich tabulek umrtnostni_tabulky <- read.csv2("D:/drctr/diplomka/souborz/umrtnostni_tabulky.csv") #install.packages("matrixStats") #library(matrixStats) #install.packages("poisson") library(poisson) library(fBasics) library(scales) ######################################################################### ######################Neasted############################################## vek_pojisteneho <- 40 # vstupni vek omega <- 106 # maximum doziti T = omega - vek_pojisteneho #T<-2 anuita <- 100000 # vyse castky pri preziti daneho roku # Model uroku level.int.rate <- 0.04 level.inf.rate <- 0.03 sigma.int.rate <- 0.008 sigma.inf.rate <- 0.008 speed.int.rate <- 0.2 speed.inf.rate <- 0.2 scenarios.1Y <- 5000 scenarios.inner <-5000 kor.coeff <- 0.8 kor <- matrix(c(1,kor.coeff,kor.coeff,1),nrow = 2) chole <- t(chol(kor)) # korelace inflace a uroku, jsou pozitivne korelovane, jak moc tot otazka # Je to spravne, ze jsou urok a inflace korelovany? # http://www.gold-eagle.com/article/gold-consequences-gibson’s-paradox umrtnost <- umrtnostni_tabulky[vek_pojisteneho+2:T, 4] #prvni rok, vim, jestli prezije umrtnost <-c(0,umrtnost[!is.na(umrtnost)]) #Prvni rok vim, jestli prezije # lamba zhorseni zdrav stavu lambda <- 25 ######################################################################## ## korelace mezi inflaci a urokem resena korelaci nahodne slozky v Vasickove modelu ##1Y scenarios #pocatecni hodnoty Inflace a uroku initial.int.rate <- 0.04 initial.inf.rate <- 0.03 #set.seed(123) rand.1Y <- matrix(data = rnorm(scenarios.1Y*2),ncol = scenarios.1Y,nrow = 2) kor <- matrix(c(1,kor.coeff,kor.coeff,1),nrow = 2) kor.hodnoty <- t(t(chol(kor))%*% rand.1Y) int.rate.1Y <- c(rep(0, scenarios.1Y)) inf.rate.1Y <- c(rep(0, scenarios.1Y)) int.rate.1Y.ant <- c(rep(0, scenarios.1Y)) inf.rate.1Y.ant <- c(rep(0, scenarios.1Y)) Rezerva.kumulativne <- rep(NA,1) for (i in 1:scenarios.1Y) { int.rate.1Y [i] <-initial.int.rate +speed.int.rate*(level.int.rate-initial.int.rate)*dt + sigma.int.rate*sqrt(initial.int.rate*dt)*kor.hodnoty[i,1] inf.rate.1Y [i] <-initial.inf.rate +speed.inf.rate*(level.inf.rate-initial.inf.rate)*dt + sigma.inf.rate*sqrt(initial.inf.rate*dt)*kor.hodnoty[i,2] } ######################################################################## #####Rezerva s Antithetic variates ######################################################################## ## Pro kazdou 1Y. hodnotu, vlastni inner scenar ## rezerva na konci Rezerva <- rep(NA,scenarios.1Y *scenarios.inner) Rezerva.ant <- rep(NA,scenarios.1Y *scenarios.inner) int.rate.2Yplus <- numeric(T-1) inf.rate.2Yplus <- numeric(T-1) int.rate.2Yplus.ant <- numeric(T-1) inf.rate.2Yplus.ant <- numeric(T-1) int.rate.2Yplus1Y <- rep(int.rate.1Y, each =scenarios.inner) inf.rate.2Yplus1Y <- rep(inf.rate.1Y, each =scenarios.inner) zdrav.stav.revize.ant<-rep(0,T) zdrav.stav.ant<-rep(0,T) ######################################################################## zacatek <-Sys.time() for (i in 1:(scenarios.1Y *scenarios.inner)) { # smycka pro kazdou jednu simulaci/scenarios.1Y *scenarios.inner # nahodna slozka pro Vasicka #set.seed(123+i) rand.2Yplus <- matrix(data = rnorm((T-1)*2),ncol = T-1 ,nrow = 2) kor.hodnoty <- t(chole%*% rand.2Yplus) # pocatecni hodnota urok/inflace int.rate.2Yplus[1] <- int.rate.2Yplus1Y[i] inf.rate.2Yplus[1] <- inf.rate.2Yplus1Y[i] int.rate.2Yplus.ant[1] <- int.rate.2Yplus1Y[i] inf.rate.2Yplus.ant[1] <- inf.rate.2Yplus1Y[i] # hodnoty uroku a inflace for (j in 2:T) { int.rate.2Yplus [j] <-int.rate.2Yplus[j-1]+speed.int.rate*(level.int.rate-int.rate.2Yplus[j-1])*dt + sigma.int.rate*sqrt(int.rate.2Yplus[j-1]*dt)*kor.hodnoty[j-1,1] inf.rate.2Yplus [j] <-inf.rate.2Yplus[j-1]+speed.inf.rate*(level.inf.rate-inf.rate.2Yplus[j-1])*dt + sigma.inf.rate*sqrt(inf.rate.2Yplus[j-1]*dt)*kor.hodnoty[j-1,2] } for (j in 2:T) { int.rate.2Yplus.ant [j] <-int.rate.2Yplus.ant[j-1]+speed.int.rate*(level.int.rate-int.rate.2Yplus.ant[j-1])*dt + sigma.int.rate*sqrt(int.rate.2Yplus.ant[j-1]*dt)*(kor.hodnoty[j-1,1]*-1) inf.rate.2Yplus.ant [j] <-inf.rate.2Yplus.ant[j-1]+speed.inf.rate*(level.inf.rate-inf.rate.2Yplus.ant[j-1])*dt + sigma.inf.rate*sqrt(inf.rate.2Yplus.ant[j-1]*dt)*(kor.hodnoty[j-1,2]*-1) } # uprava koeficienty int <- int.rate.2Yplus +1 inf <- inf.rate.2Yplus +1 int.ant <- int.rate.2Yplus.ant +1 inf.ant <- inf.rate.2Yplus.ant +1 # diskontovani t <- seq(from= 1, to = T, by = 1) int.index <- int ^ t inf.index <- inf ^ t int.index.ant <- int.ant ^ t inf.index.ant <- inf.ant ^ t #Revize nastava nahodne a vyse revize je taky nahodna #set.seed(123+i) zdrav.stav.revize<-rep(1,T) zdrav.stav.revize.ant<-rep(1,T) vyse_zmeny <- rnorm(T-1,mean=0,sd=0.3) vyse_zmeny.ant <- vyse_zmeny*-1 zdrav.stav <- floor(cumsum(rexp(T-1,lambda)))+1 for (k in 2:(T-1)) { zdrav.stav.revize[k+1] <- ifelse((zdrav.stav[k] == zdrav.stav[k-1]), zdrav.stav.revize[k] , exp(vyse_zmeny[k])) } for (k in 2:(T-1)) { zdrav.stav.revize.ant[k+1] <- ifelse((zdrav.stav[k] == zdrav.stav[k-1]), zdrav.stav.revize.ant[k] , exp(vyse_zmeny.ant[k])) } # Rezerva Rez <- sum((inf.index / int.index) * zdrav.stav.revize * (1-umrtnost))*anuita Rezerva [i] <- sum(Rez) Rez.ant <- sum((inf.index.ant / int.index.ant) * zdrav.stav.revize.ant * (1-umrtnost))*anuita Rezerva.ant [i] <- sum(Rez.ant) #bez.na <- subset(Rezerva, !is.na(Rezerva)) #Rezerva.kumulativne[i]<-mean(bez.na) } Rezerva_tot <- c(Rezerva.ant,Rezerva) Rezerva_tot <- as.vector(rbind(Rezerva.ant,Rezerva)) ######################################################################## ########### konec <-Sys.time() print(konec - zacatek) zaloha500_5000<-Rezerva_tot prumery_neasted <-running(rezerva_zaloha, width=scenarios.1Y, by=scenarios.1Y) ##alternativne sekvence <- seq(scenarios.inner,scenarios.inner*scenarios.1Y,by=scenarios.inner) prumery <- rep(NA,scenarios.1Y) pomocna <- 1 i <- 5000 k<-0 for (i in sekvence) { prumery[pomocna] <- mean(Rezerva.zaloha2[k:i]) k<-i+1 pomocna<- pomocna+1 } sekvence[1] mean(Rezerva.zaloha2[5001:10000]) sd(prumery) skewness(prumery) kurtosis(prumery) quantile(prumery,c(0.99)) prumery[1] Rezerva.zaloha2[1] mean(Rezerva.zaloha2) ############################### ############################### ###Graf CIR x <- seq(1,66,1) plot(x,int.rate.2Yplus,lwd = 2, yaxt="n",type="l",col="red",xlab="t", ylab="Úrok/Inflace", ylim=c(0.02, 0.05)) lines(x,inf.rate.2Yplus,col="blue",lwd = 2 ) legend(50,0.05, # umisteni legendy v grafu c("Úrok","Inflace"), # text v legende lty=c(1,1), # symboly do legendy lwd=c(2,2),col=c("red","blue"),cex = 1) # barva a sirka krivek returns <- runif(100) yticks_val <- seq(0.02,0.1,0.01) axis(2, at=yticks_val, lab=percent(yticks_val)) ############################### ############################### ###Graf Poisson zdrav.stav.revize<-rep(0,T) zdrav.stav <- floor(cumsum(rexp(T-1,lambda)))+1 for (k in 2:(T-1)) { zdrav.stav.revize[k+1] <- ifelse((zdrav.stav[k] == zdrav.stav[k-1]), 0 , rnorm(1,sd=0.3)) } zdrav.stav.revize <- cumsum(zdrav.stav.revize)+1 plot(zdrav.stav.revize,type= "l",xlab="t",ylab="Zdravotní stav",yaxt="n",ylim=c(0,2)) yticks_val <- seq(0,2,1) axis(2, at=yticks_val, lab=percent(yticks_val)) ############################### ############################### ###Histogram rezervy hist(prumery/1000000, breaks = 20 ,xlab="Nejlepší odhad za 1 rok (v mil)", ylab="Četnost",main = NULL) yticks_val <- seq(3,9,1) axis(1, at=yticks_val, lab=(yticks_val)) ?axis max(Rezerva) min(Rezerva) mode(Rezerva) median(Rezerva) mean(Rezerva) ############################### ############################### ###Graf neasted monte carlo CIR <- matrix(0,ncol = 59,nrow = 5) initial <- matrix(0.04,ncol = 1,nrow = 5) CIR <- cbind(initial,CIR) CIR [,2] <-0.04+ sigma.int.rate*sqrt(0.04*dt)*rnorm(1) for (i in 1:5) { for (j in 3:60) { CIR [i,j] <-CIR[i,j-1]+speed.int.rate*(level.int.rate-CIR[i,j-1])*dt + sigma.int.rate*sqrt(CIR[i,j-1]*dt)*rnorm(1) } } CIR proces<- function(){ CIR <- rep(0,66) CIR[1] <- (0.04) for (j in 2:61) { CIR [j] <-CIR[j-1]+speed.int.rate*(level.int.rate-CIR[j-1])*dt + sigma.int.rate*sqrt(CIR[j-1]*dt)*rnorm(1) } return(CIR) } time <- rep(0,100) for (i in 1:100) { time[i] <- system.time(proces())[3] } mean(time) plot(CIR[1,],lwd = 2, yaxt="n",type="l",col="red",xlab="t", ylab="Úrok", ylim=c(0.032, 0.046)) lines(CIR[2,],col="blue",lwd = 2 ) lines(CIR[4,],col="green4",lwd = 2 ) lines(CIR[5,],col="yellow2",lwd = 2 ) lines(CIR[3,],col="black",lwd = 2 ) yticks_val <- seq(0.03,0.05,0.005) axis(2, at=yticks_val, lab=percent(yticks_val)) ?axis ############################################### ############################################### ### Vypocetni cas 1 smycky proces<- function(){ CIR <- rep(0,66) CIR[1] <- (0.04) for (j in 2:60) { CIR [j] <-CIR[j-1]+speed.int.rate*(level.int.rate-CIR[j-1])*dt + sigma.int.rate*sqrt(CIR[j-1]*dt)*rnorm(1) } return(CIR) } time <- rep(0,1000) for (i in 1:1000) { t1 <- Sys.time() time[i] <- proces() t2 <- Sys.time() time[i] <- t2-t1 } mean(time) warnings()