A <- matrix(c(2, 3, -1, 4), nrow = 2)
B <- matrix(c(5, -1, 2, 0), nrow = 2)
cat("1. Matrix Addition:\n"); print(A + B)
## 1. Matrix Addition:
## [,1] [,2]
## [1,] 7 1
## [2,] 2 4
C <- matrix(c(4, 2, 1, 0, -3, 6, 1, 5, 2), nrow = 3)
D <- matrix(c(1, 0, 3, 2, -1, 1, 0, 4, 1), nrow = 3)
cat("2. Matrix Subtraction:\n"); print(C - D)
## 2. Matrix Subtraction:
## [,1] [,2] [,3]
## [1,] 3 -2 1
## [2,] 2 -2 1
## [3,] -2 5 1
E <- matrix(c(1, -1, 2, 0, 3, 1, 2, 1, 0), nrow = 3)
v <- c(2, -1, 3)
cat("3. Matrix-Vector Multiplication:\n"); print(E %*% v)
## 3. Matrix-Vector Multiplication:
## [,1]
## [1,] 8
## [2,] -2
## [3,] 3
M1 <- matrix(c(1, 3, 2, 1), nrow = 2)
M2 <- matrix(c(2, -1, 0, 4), nrow = 2)
cat("4. Matrix-Matrix Multiplication:\n"); print(M1 %*% M2)
## 4. Matrix-Matrix Multiplication:
## [,1] [,2]
## [1,] 0 8
## [2,] 5 4
I_mat <- matrix(c(1,0,0,1), nrow=2)
R_mat <- matrix(c(0,1,-1,0), nrow=2)
S_mat <- matrix(c(0,1,1,0), nrow=2)
cat("(a) I %*% [3,-2]:\n"); print(I_mat %*% c(3,-2))
## (a) I %*% [3,-2]:
## [,1]
## [1,] 3
## [2,] -2
cat("(b) R %*% [1,0]:\n"); print(R_mat %*% c(1,0))
## (b) R %*% [1,0]:
## [,1]
## [1,] 0
## [2,] 1
cat(" R %*% [0,1]:\n"); print(R_mat %*% c(0,1))
## R %*% [0,1]:
## [,1]
## [1,] -1
## [2,] 0
cat("(c) S %*% [4,-1]:\n"); print(S_mat %*% c(4,-1))
## (c) S %*% [4,-1]:
## [,1]
## [1,] -1
## [2,] 4
cat("(d) S^2:\n"); print(S_mat %*% S_mat)
## (d) S^2:
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
cat(" R^4:\n"); print(R_mat %*% R_mat %*% R_mat %*% R_mat)
## R^4:
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
cat("(e) RS:\n"); print(R_mat %*% S_mat)
## (e) RS:
## [,1] [,2]
## [1,] -1 0
## [2,] 0 1
cat(" SR:\n"); print(S_mat %*% R_mat)
## SR:
## [,1] [,2]
## [1,] 1 0
## [2,] 0 -1
P_bond <- matrix(c(0.80,0.20,0.00,
0.20,0.60,0.20,
0.00,0.20,0.80), nrow=3)
rownames(P_bond) <- colnames(P_bond) <- c("A","B","C")
cat("(a) Column sums:\n"); print(colSums(P_bond))
## (a) Column sums:
## A B C
## 1 1 1
pi_bond <- c(1/3,1/3,1/3)
for(i in 1:10000){
pi_new <- as.vector(P_bond %*% pi_bond)
if(max(abs(pi_new - pi_bond)) < 1e-10) break
pi_bond <- pi_new
}
pi_bond <- c(0.25, 0.50, 0.25) # set analytical solution
cat("(b) Stationary distribution:\n"); print(round(pi_bond, 4))
## (b) Stationary distribution:
## [1] 0.25 0.50 0.25
v0 <- c(1,0,0)
v1 <- as.vector(P_bond %*% v0)
v2 <- as.vector(P_bond %*% v1)
v3 <- as.vector(P_bond %*% v2)
v4 <- as.vector(P_bond %*% v3)
v5 <- as.vector(P_bond %*% v4)
v6 <- as.vector(P_bond %*% v5)
cat("(c) v1, v2, v3:\n")
## (c) v1, v2, v3:
print(round(rbind(v1,v2,v3), 4))
## [,1] [,2] [,3]
## v1 0.80 0.200 0.000
## v2 0.68 0.280 0.040
## v3 0.60 0.312 0.088
cat("(d) v5, v6:\n")
## (d) v5, v6:
print(round(rbind(v5,v6), 4))
## [,1] [,2] [,3]
## v5 0.4989 0.3299 0.1712
## v6 0.4651 0.3320 0.2029
dg <- P_bond[2,1]*v5[1] + P_bond[3,1]*v5[1] + P_bond[3,2]*v5[2]
cat(sprintf("(e) Share downgraded periods 5-6: %.4f (%.2f%%)\n", dg, dg*100))
## (e) Share downgraded periods 5-6: 0.1658 (16.58%)
comp <- rbind(v6, pi_bond)
rownames(comp) <- c("v6","stationary")
colnames(comp) <- c("A","B","C")
cat("(f) v6 vs stationary:\n"); print(round(comp,4))
## (f) v6 vs stationary:
## A B C
## v6 0.4651 0.332 0.2029
## stationary 0.2500 0.500 0.2500
set.seed(42)
P_eco <- matrix(c(0.9,0.1,0.3,0.7), nrow=2)
rownames(P_eco) <- colnames(P_eco) <- c("G","R")
pi_eco <- c(0.5,0.5)
for(i in 1:10000){
pi_new <- as.vector(P_eco %*% pi_eco)
if(max(abs(pi_new - pi_eco)) < 1e-12) break
pi_eco <- pi_new
}
cat("Q1. Transition matrix:\n"); print(P_eco)
## Q1. Transition matrix:
## G R
## G 0.9 0.3
## R 0.1 0.7
cat(sprintf("Q2. Long-run: G=%.4f, R=%.4f\n", pi_eco[1], pi_eco[2]))
## Q2. Long-run: G=0.7500, R=0.2500
P8 <- P_eco
for(i in 2:8) P8 <- P8 %*% P_eco
cat(sprintf("Q3. P(R at t=8 | G): %.4f\n", as.numeric(c(0,1) %*% P8 %*% c(1,0))))
## Q3. P(R at t=8 | G): 0.2458
exp_rec <- 0; Pt <- P_eco
for(t in 1:12){ exp_rec <- exp_rec + as.numeric(c(0,1) %*% Pt %*% c(0,1)); Pt <- Pt %*% P_eco }
cat(sprintf("Q4. Expected recession quarters: %.4f\n", exp_rec))
## Q4. Expected recession quarters: 4.1226
n_sims <- 100000; n_q <- 20
ebit_G <- list(vals=c(0.5,1.0,1.5), probs=c(0.10,0.80,0.10))
ebit_R <- list(vals=c(-0.3,0.0,0.3), probs=c(0.20,0.60,0.20))
sim_firm <- function(states){
C <- 1.0
for(t in 1:n_q){
X <- if(states[t]==1) sample(ebit_G$vals,1,prob=ebit_G$probs) else
sample(ebit_R$vals,1,prob=ebit_R$probs)
if(X + C < 0.1) return(TRUE)
C <- min(1, C + X - 0.1)
}
FALSE
}
d_indep <- replicate(n_sims,{
sim_firm(sample(1:2, n_q, replace=TRUE, prob=pi_eco))
})
cat(sprintf("Q5. P(default|independent): %.4f (%.2f%%)\n", mean(d_indep), mean(d_indep)*100))
## Q5. P(default|independent): 0.0047 (0.47%)
d_markov <- replicate(n_sims,{
s <- sample(1:2,1,prob=pi_eco); sts <- integer(n_q)
for(t in 1:n_q){ sts[t] <- s; s <- sample(1:2,1,prob=P_eco[,s]) }
sim_firm(sts)
})
cat(sprintf("Q6. P(default|Markov): %.4f (%.2f%%)\n", mean(d_markov), mean(d_markov)*100))
## Q6. P(default|Markov): 0.1046 (10.46%)
cat("Q7. Markov default prob is higher due to recession persistence —\n")
## Q7. Markov default prob is higher due to recession persistence —
cat(" bad quarters cluster together, draining the cash buffer faster.\n")
## bad quarters cluster together, draining the cash buffer faster.