# -------------------------------------------------------------------------------------------------------------------------------- # FCI indicator: Plasil, Seidler, Konecny, Hlavac - In the quest of measuring the financial cycle, WP CNB # -------------------------------------------------------------------------------------------------------------------------------- # The code comes with absolutely no warranty. Adjustments to the lines 11, 12 and 62 must be done in order to reflect country-specific features. # ------------------------------------------------------------------------------------------------------------------------------------ # The code makes use of the ks package (kdce function): Make sure that it is installed library(ks) # Data load and manipulation --------------------------------------------------------------------------------------------------- RawData <- data.matrix(read.csv2("input.csv")) # Data load (variables should be stored in columns with variable names in the first row) Multiply.columns <- c(1,1,1,1,1,-1,-1,1,-1) # premultiply some variables by -1 (if necessary) so that ... RawData <- RawData %*% diag(Multiply.columns) # ...values around 0 represent the trough and values around 1 represent the peak of the cycle nVar <- ncol(RawData) nObs <- nrow(RawData) Indicators <- apply(RawData,2,function(x) kcde(x, eval.points = x)$estimate) # map the variables onto 0-1 interval using kernel estimate of the distribution function # Time-varying correlations --------------------------------------------------------------------------------------------------- lambda <- 0.94 # EWMA parameter # correlations at time t=1 ewma.weights <- (lambda^(0:(nObs-1))) / sum(lambda^(0:(nObs-1))) CovInit <- cov.wt(Indicators, wt= ewma.weights)$cov # pairwise time-varying correlations CovArray <- array(NA, c(nVar,nVar,nObs)) CorArray <- array(NA, c(nVar,nVar,nObs)) CovArray[,,1] <- CovInit for (i in 1:nVar) { for (j in 1:nVar) { for (q in 2:nObs) { CovArray[i,j,q] <- lambda*CovArray[i,j,(q-1)] + (1-lambda)*(Indicators[q,i]-0.5)*(Indicators[q,j]-0.5) # EWMA estimate of the corr. matrix } } } for (q in 1:nObs) { CorArray[,,q] <- diag(1/sqrt(diag(CovArray[,,q]))) %*% CovArray[,,q] %*% diag(1/sqrt(diag(CovArray[,,q]))) } # FCI estimate (Sick was the original name for the FCI) ---------------------------------------------------------------------------------------- Sick.values <- numeric() Sick.theor <- numeric() Sick.decomp <- matrix(NA, nObs,nVar) Indi.weights <- c(0.35, 0.27, 0.09, 0.08, 0.07, 0.05, 0.05, 0.02, 0.02) # values of the fixed weights w for (q in 1:nObs) { Sick.values[q] <- (Indicators[q,]*Indi.weights) %*% CorArray[,,q] %*% (Indicators[q,]*Indi.weights) # actual values of the FCI Sick.theor[q] <- (Indicators[q,]*Indi.weights) %*% matrix(1,nVar,nVar) %*% (Indicators[q,]*Indi.weights) # theoretical values under the assumption of perfect correlations Sick.decomp[q,] <- (Indicators[q,]*Indi.weights) %*% matrix(1,nVar,nVar) * (Indicators[q,]*Indi.weights) # contributions of individual variables to the FCI } plot.ts(Sick.values) # Chart # Appendix #------------------------------------------------------------------------------------------------------------------------------- # Simulation to find out/calibrate optimal weights w #------------------------------------------------------------------------------------------------------------------------------- y <- CreditLoss # insert a dependent variable of interest, e.g.(future) credit losses (at time t+6 quarters) nSim <- 30000 # number of simulations rmse.res <- matrix(NA, nSim, (nVar+1)) Sick.values <- numeric() Sick.theor <- numeric() for (i in 1:nSim) { Indi.weights <- sort(runif(nVar),decreasing = T) Indi.weights <- Indi.weights/sum(Indi.weights) for (q in 1:nObs) { Sick.values[q] <- (Indicators[q,]*Indi.weights) %*% CorArray[,,q] %*% (Indicators[q,]*Indi.weights) } x <- Sick.values # use the FCI with simulated weights in a "predictive" regression (you may want to use a more complicated model) rmse.res[i,] <- c(Indi.weights, sqrt((t(y- lm(y~x)$fitted.values) %*% (y- lm(y~x)$fitted.values))/nObs)) # RMSE (other criteria are possible) } rmse.res <- rmse.res[order(rmse.res[,10]),] # show 10 best performing models with respect to the RMSE of the model # ------------------------------------------------------------------------------------------------------------------------------