R - Monte Carlo simulation for Credit Risk (PD) Stress Testing Scenarios

################################################################################
# Stress Test Forecasting Model
# Monte Carlo Simulations
#
# Inputs:
# 1. Base portfolio, customers list with baseline probability of default.
# 2. Transition probability matrices, 1 per year, per scenario.
#
# Output:
# 1. Forecasted list of obligors defaulted per year/scenario... with an
#    expected default frequency.
#
#
################################################################################




####
#### 1. Data import
####

# 1. Set default working path
  setwd('S:/TMP Directory')

# 2. Import all obligors
# Columns required: Customer Id, Base Probability of Default
obl = read.csv('01 - Customers with PD.csv',header=T)

# 3. ORs
# We might move this out into an external file later on.
ORs = c(0.0001, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.0012, 0.30, 1)

# 4. TPMs
# Use the same oTPM for the three years for base case.
# Use a different TPM for each year when calculating other scenarios such as Stress1 and Stress2.
oTPM1 = read.csv( '02.01 - TPM USA Stress 2 2010.csv', header = T )
oTPM2 = read.csv( '02.02 - TPM USA Stress 2 2011.csv', header = T )
oTPM3 = read.csv( '02.03 - TPM USA Stress 2 2012.csv', header = T )




####
#### 2. declarations and initializations
####


# 2.1. Define the Outupt Data Structure
output = data.frame( run=0, yr=0, rowid=0, customerId='a', basePD=0 )

# 2.2. Number of montecarlo simulations
numsims = 50

# 2.3. Number of years.
# Keep in mind that you will need to have loaded and to apply a different TPM
# for each of the years.
numyrs = 3




####
#### 3. Run the simulations
####
print( '3. Simulation...' )

# For each simulation
for( run in 1:numsims ){

# For each year
for( yr in 1:numyrs ){

# Apply the first TPM matrix to the first year
if ( yr == 1 ) {
strt = as.vector( obl$basePD )
tmp = strt
oTPM_local = oTPM1
} #yr1

# Apply the second TPM matrix to the second year
if( yr == 2 ) {
strt = tmp
oTPM_local = oTPM2
} #yr2

# Apply the third TPM to the third year
if( yr == 3 ) {
strt = tmp
oTPM_local = oTPM3
} #yr3

# Random number genarator draws
unif_draws = runif( length( tmp ) )
for( i in 1:length( ORs ) ){
breaks = oTPM_local[ oTPM_local$startPD == ORs[i] & oTPM_local$Prob != 0, 'LowerBound' ]
ORLevels = oTPM_local[ oTPM_local$startPD == ORs[i] & oTPM_local$Prob != 0, 'UpperBound' ]

#Just to be safe, let's sort them
breaks = breaks[ order( ORLevels ) ]
ORLevels = ORLevels[ order( ORLevels ) ]
tmp[ strt == ORs[i] ] = as.numeric( as.character( cut( x = unif_draws[strt == ORs[i]], c(breaks,1), ORLevels ) ) )
} #draw

# Only output if this transition to ORD (if it defaulted)
if( sum( strt != 1 & tmp == 1 ) > 0 ){
output = rbind(   output, data.frame( run=run, yr=yr,
rowid=obl[strt!=1 & tmp==1,'rowid'],
customerId=obl[strt!=1 & tmp==1,'customerId'],
basePD=obl[strt!=1 & tmp==1,'basePD'])   )
} #if defaulted

} #each year

} #each run/sim




####
#### 4. Output
####
print( '4. About to output... ' )
output = output[ 2 : nrow(output) , ]
table( output[ , c( 'run', 'yr' ) ] )
write.csv( output,'Defaulters.csv' )