2. Long run proportion of time in each rating
To determine the long run proportion of time within each rating, we note that the net transition rate into one category is equal to the net transition rate out of the same category. This can be translated into the following equations:
and
where are the transition rates from rating k to rating j, is the transition rate out of rating j, while is the long run proportion of time in rating k. The values of and for every rating can be found in the estimated transition rate matrix. The second equation is from noting that the proportion of time in each rating should add up to 1. These equations can be solved simultaneous using matrix operations, and MATLAB is once again used to obtain the results. The command lines and code used to run this procedure can be found in Appendix 1.2.
The long run proportion of time in each rating is shown in the below table:
These proportions can be verified by calculating the probability transition matrix P(t) for large values of t. It should be noted that in the long run, sovereigns will most likely spend a significant proportion of their time in AAA rating (at around 76% of time). This is because the transition rate from AAA rating to other ratings is much lower than for other ratings.
3. Non constant transition rates
Now consider the scenario where the transition rates start increasing by 10% after 1999. This can be represented by multiplying the transition rate matrix R by 1.1 for every year it increases. We are asked to calculate the probabilities of a sovereign from each category being rated as default in 3, 5 and 10 years as at the beginning of 1995.
Since the transition matrix will only increase after 1999, 5 years from 1995, the 3 year and 5 year probabilities can be calculated same as before, by calculating P(t) based on . MATLAB is used to find the exponential of the matrix R and the appropriate probabilities are extracted from the matrices P(3) and P(5).
For the 10 year probabilities, the solution of P(t) will change every year after 5 years due to changes in transition rates matrix. For example, the 6th year (1999-2000), the solution is now, where t is time after 1999 and before 2000. The Chapman Kolmogorov equations, which areare utilised to account for this problem, by calculating P(5) based on the unchanged transition matrix for years between 1995 and 1999, while calculating P(1) for each year where the transition rates matrix R increases for the 5 years after 1999, so that P(10) is just the product of all these matrices. Once again, MATLAB is used to calculate the exponential of matrices.
The table below shows the probabilities of a sovereign from each category being rated as default (SD) in 3, 5 and 10 years as at the beginning of 1995. The code and further results can be found in the Appendix 1.3.
It can be observed that the probability of a sovereign initially rated BB or above being rated default increases over time, due to greater time frame for credit events and hence rating migration to occur. Also the increase in transition rate matrix means that there will be more transitions between credit ratings, meaning sovereigns starting from a rating close to default are more likely to move upwards than usual.
Part II1. a) Consider the stock price modelled by the geometric Brownian motion process:
For each of the 2 methods suggested, we simulated S(t) for t=150000 with increments of 1/12 for time and then calculated ln(S(t))/t for each value of t.
It can be noted from the above results, both methods will cause the function ln(S(t))/t will converge to the limit which is approximately 0.0558077, which is close to the true value of p of 0.055. The true value can be found based on the property of the probability distribution of ln(S(t)).
We note that: as S(0)=1, , and
Using the Delta Method, we can deduce the distribution of ln(S(t))/t as:
Taking limits as t goes towards infinity, we note that:
So as t goes to infinity, the variance will reduce to 0, meaning ln(S(t))/t will approach the mean which is 0.055. Thus the true value of p is 0.055 and verified the results of the simulation.
b) To find the value of t to obtain 2 digit accuracy with 95% probability, we find t such that:
This is equivalent to finding:
Rearranging gives
Noting that from part a) that , where p=v, we can deduce that . This implies that , and noting that and , we deduce that .
So t must be around 4571822 to obtain two-digit accuracy of the true value of p with approximately 95% probability.
c) Based on the values of S(t) simulated from both methods, the value of can also be simulated using both methods. The code can be found in Appendix 2.1.3.
It can be observed from the above graphs that the function does not converge to a limit. The function will continually fluctuate randomly, following a distribution similar to Gamma random distribution.
The behaviour can be verified by analysing the properties of probability distribution of , which we obtain from part (a) as:
This implies that. It can then be shown that
where denotes the Gamma distribution with parameters 0.5 and 0.18. The full derivation can be found in the Appendix 2.1.2. Clearly the behaviour of is independent of t, and thus the function does not approach a limit as t tends to infinity.
2. Weekly exchange rate analysis – US Dollar and pound sterling (1980-1988)
Before constructing a stochastic model to fit the data, the data is plotted as a times series, and the sample autocorrelation function (ACF) and sample partial autocorrelation function (PACF) are plotted to see if there is any observable trend or seasonality.
It can be noted that from the above plot of the time series of weekly exchange rates, there does not seem to be any significant trend or seasonality. From the sample ACF plot, the correlations decrease rather slowly and steadily from 1, indicating the exchange rates are not stationary. Thus differencing the data may be advisable before fitting a model. The time series of the differenced data is then plotted, as well as the sample ACF and sample PACF plots.
The above plots indicate the differenced data for exchange rates are stationary, as the ACF rapidly goes towards zero after lag 1, while the partial ACF is within the , with This suggests an ARIMA(0,1,0) model may be an appropriate model for exchange rates.
When the ARIMA(0,1,0) was fitted, there are certain criteria that needs to be analysed. The AIC (Akaike Information Criteria) for this model is -2035.93 (Other information can be found in the Appendix). To see that this model minimises the AIC, the AIC from the ARIMA(0,1,1) and ARIMA(1,1,0) models are found and compared with the AIC for the proposed model. The AIC from ARIMA(0,1,1) model is -2035.74 while the AIC from ARIMA(1,1,0) model is -2035.79. The AIC of the proposed ARIMA(0,1,0) model is clearly smaller than the AIC of the other 2 models.
Another way to assess the model is to assess the behaviour of the residuals from the fitted model. The standardised residuals, sample ACF of residuals and p values of Ljung-Box (Portmanteau) statistics are plotted on the graph below.
From the residual plot, it can be observed that the residuals form a good approximation to a white noise process, with no real seasonal pattern in fluctuation or any trend in the residuals. The ACF of Residuals plot shows that most of the sample autocorrelations lie within the 95% confidence interval indicated by the dotted lines. The p values for the Ljung-Box (Portmanteau) statistics are all outside the rejection region, which means that the null hypothesis, which states the residuals are white noise with zero mean and variance 1, is accepted. The three above plots all indicate that ARIMA(0,1,0) is an appropriate model of the time series of exchange rates.
Based on the analysis of the residuals for the fitted model and the AIC, it can be concluded that the exchange rates can be fitted using an ARIMA(0,1,0) model, which implies , where is the exchange rate at time t, B is the backshift operator and are uncorrelated white noise. This is equivalent to, which suggests that exchange rates follow a random walk process, which is typical for most financial data. Thus, most economists will try to simulate exchange rate movements using the random walk (ARIMA(0,1,0)) process.
This analysis was performed in R. The code and further results can be found in the Appendix.
Appendix
1. Part I
1.1 – Code for Question 1
First, the matrix R was entered into MATLAB and stored as “R”. Then P(t) is calculated for 100 equally spaced points between time 0 and time 10. The relevant entries are collected and stored and then plotted.
Create vectors to store transition probabilities for each rating
p=linspace(1,10,100); For AAA
q=linspace(1,10,100); For AA
r=linspace(1,10,100); For A
s=linspace(1,10,100); For BBB
format long
for t=1:100
A=expm(R*(t/10)); Function to find exponential of matrices, noting P(0) is identity matrix
Collects entries which are relevant to each rating and store them
p(t)=A(73);
q(t)=A(74);
r(t)=A(75);
s(t)=A(76);
end
t=linspace(1,100);
Plot probabilities of default over time for all 4 ratings
figure(1)
plot(t/10,p,'-r',t/10,q,'-b',t/10,r,'-k',t/10,s,'-g')
t/10 scales time to between 0 and 10
xlabel('Time (years from 1995)')
ylabel('Probability of default rating')
title('Probability of default for AAA, AA, A, BBB')
grid on
Plot probabilities of default over time for AAA, AA and A
figure(2)
plot(t/10,p,'-r',t/10,q,'-b',t/10,r,'-k')
t/10 scales time to between 0 and 10
xlabel('Time (years from 1995)')
ylabel('Probability of default rating')
title('Probability of default for AAA, AA, A')
grid on
Plot probabilities of default over time for AAA, AA
figure(3)
plot(t/10,p,'-r',t/10,q,'-b')
t/10 scales time to between 0 and 10
xlabel('Time (years from 1995)')
ylabel('Probability of default rating')
title('Probability of default for AAA, AA')
grid on
Plot probabilities of default over time for AAA
figure(4)
plot(t/10,p,'-r')
t/10 scales time to between 0 and 10
xlabel('Time (years from 1995)')
ylabel('Probability of default rating')
title('Probability of default for AAA')
grid on
In order to find the long run probabilities:
Enter matrix to solve by deleting last row of R, transposing it and then adding a row vectors of 1 to add the condition that all probabilities will add to one.
m =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
-0.0059 0.0544 0 0 0 0 0 0 0
0.0059 -0.0816 0.0251 0 0 0 0 0 0
0 0.0272 -0.0564 0.0945 0 0 0 0 0
0 0 0.0313 -0.1323 0.0691 0 0 0 0
0 0 0 0.0315 -0.1487 0.1022 0 0 0
0 0 0 0.0063 0.0797 -0.1981 0.3019 0 0.9677
0 0 0 0 0 0.0767 -0.9811 0 0.1935
0 0 0 0 0 0 0.3774 -3.1579 0
Create column vector so that mP=b where P is the long run probability vector
>> b=[1;0;0;0;0;0;0;0;0]
>> m\b Solve for P in mP=b
ans =
0.7642
0.0829
0.0898
0.0298
0.0163
0.0145
0.0013
0.0002
0.0010
2.1 Part II Question 12.1.1
Simulates values of S(t) for 150000 values of t, noting each value of t increments by 1/12.
t=1
Generate vectors used to store values of S(t) for both methods
S1=seq(1:150000)
S1[1]=1 Noting S(0)=1
S2=seq(1:150000)
S2[1]=1 Noting S(0)=1
Generate vectors used to store values of lnS(t)/t for both methods
Q1=seq(1:150000)
Q2=seq(1:150000)
Generate vectors used to store values of [(lnS(t)-pt)^2]/t for both methods
P1=seq(1:150000)
P2=seq(1:150000)
while (t<150001) { Loops until 150000 values are generated
u=0.1
o=0.3
dt=1/12
v=u-o^2/2
e1=rnorm(1,0,1)
if (t<150000) {
S1[t+1]=(1+u*dt+o*e1*sqrt(dt))*S1[t]
S2[t+1]=exp(v*dt+o*e1*sqrt(dt))*S2[t]
}
Q1[t]=12/t*log(S1[t]) Need to multiply by 12 as time increments as 1/12
Q2[t]=12/t*log(S2[t])
P1[t]=12/t*(log(S1[t])-0.055*t/12)^2
P2[t]=12/t*(log(S2[t])-0.055*t/12)^2
t=t+1
}
Plots function for 1(a) Method 1
plot(Q1, type="l", lwd=3, xlab="Time", ylab="(ln[S(t)])/t", main="Result from Method 1")
Plots function for 1(a) Method 2
plot(Q2, type="l", lwd=3, xlab="Time", ylab="(ln[S(t)])/t", main="Result from Method 1")
>Q1[150000] = 0.0558077 ←Used as an estimate of the limit
Plots function for 1(c) Method 1
plot(P1, type="l", lwd=3, xlab="Time", ylab="((ln[S(t)]-p)^2)/t", main="Result from Method 1")
Plots function for 1(c) Method 2
plot(P2, type="l", lwd=3, xlab="Time", ylab="((ln[S(t)]-p)^2)/t", main="Result from Method 2")
2.1.2 – Derivation of Gamma distribution of the function
Consider and hence . So:
Let , so then:
by noting the cumulative distribution function of standard normal
Thus the probability density function of can be found by differentiating with respect to x. Using the Fundamental theorem of calculus, this gives:
where , which is in the form of distribution.
2.2 Part II Question 2Code is in Blue
Results are in Red
Comments are in Green
Collects the exchange rate data from a CSV file
exchange <- read.csv("C:/Users/Chung-Yu/Documents/Uni/Actuarial Studies/ACTL2003/Assignment/exchange.csv",header=T)
exchange<-ts(exchange[,1],start=1980,freq=365/7)
Plots time series of data as well as sample ACF and PACF
par(mfrow=c(2,2))
plot(exchange)
acf(exchange)
acf(exchange, type="partial")
Plots the differenced data
par(mfrow=c(2,2))
plot(diff(exchange))
acf(diff(exchange))
acf(diff(exchange), type="partial")
fit=arima(exchange, order=c(0,1,0)) fit ARIMA(0,1,0) model
tsdiag(fit) Shows residual plots for fitted ARIMA(0,1,0) model, used for residual analysis
> fit=arima(exchange, order=c(0,1,0))
> fit
Shows AIC for ARIMA(0,1,0) model
Call:
arima(x = exchange, order = c(0, 1, 0))
sigma^2 estimated as 0.0007593: log likelihood = 1018.97, aic = -2035.93
The next 2 ARIMA models are fitted to compare the AIC and to see that the AIC for ARIMA(0,1,0) is optimal.
> fit=arima(exchange, order=c(0,1,1))
> fit
Call:
arima(x = exchange, order = c(0, 1, 1))
Coefficients:
ma1
-0.0610
s.e. 0.0449
sigma^2 estimated as 0.0007563: log likelihood = 1019.87, aic = -2035.74
> fit=arima(exchange, order=c(1,1,0))
> fit
Call:
arima(x = exchange, order = c(1, 1, 0))
Coefficients:
ar1
-0.0629
s.e. 0.0461
sigma^2 estimated as 0.0007563: log likelihood = 1019.89, aic = -2035.79