#(c) Pavol Jurca, 2014 rm(list = setdiff(ls(), lsf.str())) Data <- read.table(file="D:/Documents/Skola/Vyuka/Risk management/Slidy/3 Mixture distributions/Exercises/fx_rates.txt",,header=TRUE, dec = ",") L=nrow(Data) Returns <- log(Data[2:L,2:6])-log(Data[1:L-1,2:6]) #MIXTURE OF 2 NORMAL DISTRIBUTION, 1 MAIN RISK FACTOR require(mclust) obj <- densityMclust(Returns[,1],G=2,control=emControl(eps=1e-20,tol=1e-6,itmax=100)) Proportions <- obj$parameters$pro Means <- obj$parameters$mean Volatilities <- sqrt(obj$parameters$variance$sigmasq) hectic <- if(Volatilities[2]>Volatilities[1]) 2 else 1 #Calculating parameters for the primary risk factor print("Calculating parameters for the primary risk factor") quiet <- 3-hectic; omega <- Proportions[hectic] mu2 <- Means[hectic] mu1 <- Means[quiet] sigma2 <- Volatilities[hectic] sigma1 <- Volatilities[quiet] #Calculating probablities that an observation is from hectic period h <- (omega*dnorm(Returns[,1],mu2,sigma2))/((1-omega)*dnorm(Returns[,1],mu1,sigma1)+omega*dnorm(Returns[,1],mu2,sigma2)) plot(Returns[,1],h,col="red",type = "p",lwd=1) #Calculating parameters for the secondary risk factor print("Calculating parameters for peripheral risk factor") scenarios=c(0.3,0,0,0,0,0) for(i in 2:6){ N <- length(Returns[,i]) mu_y2 <- sum(h*Returns[,i])/sum(h) mu_y1 <- sum((1-h)*Returns[,i])/sum(1-h) sigma_y2 <- sqrt(sum(h*(Returns[,i]-mu_y2*rep(1,N))^2)/sum(h)) sigma_y1 <- sqrt(sum((1-h)*(Returns[,i]-mu_y1*rep(1,N))^2)/sum(1-h)) rho_y2 <- sum(h*(Returns[,1]-mu2*rep(1,N))*(Returns[,i]-mu_y2*rep(1,N)))/(sigma2*sigma_y2*sum(h)) rho_y1 <- sum((1-h)*(Returns[,1]-mu1*rep(1,N))*(Returns[,i]-mu_y1*rep(1,N)))/(sigma1*sigma_y1*sum(1-h)) scenarios[i] <- rho_y2*(scenarios[1]-mu2)*sigma_y2/sigma2+mu_y2 print(c("mu_y2","mu_y1","sigma_y2","sigma_y1","rho_y2","rho_y1","scenarios[i]")) print(c(mu_y2,mu_y1,sigma_y2,sigma_y1,rho_y2,rho_y1,scenarios[i])) pause <- readline() } #TEST OF THE SPECIFICATION OF THE MODEL FOR 1 MAIN RISK FACTOR print("TEST OF THE SPECIFICATION OF THE MODEL FOR 1 MAIN RISK FACTOR") temp_obj <- densityMclust(Returns[,1],G=2,control=emControl(eps=1e-20,tol=1e-6,itmax=100)) LogL2 <- temp_obj$loglik temp_obj <- densityMclust(Returns[,1],G=1,control=emControl(eps=1e-20,tol=1e-6,itmax=100)) LogL1 <- temp_obj$loglik Likelihood_ratio <- -2*(LogL1-LogL2) Likelihood_ratio ch2_quantil <- qchisq(.95, df=3) ch2_quantil