############ Exercises from chapter 13 ########
## 18.04.18 ## 13.17-13.21##
#### 13.17
#a)
#Compute a general function which generates m sequences of a Gaussian AR(2) up to time K.
#Initial values are 0.
#m = number of realizations
#a = vector of AR-coefficients
#sigma = standard deviation of the errors
ex1317=function(m,K,sigma,a){
#start matrix with initial condition 0
X=matrix(0,K+2,m)
for (k in 1:K) {
X[k+2,]=a[1]*X[k+1,]+a[2]*X[k,]+sigma*rnorm(m)
}
list(X=X)
#X = (K+2)xm matrix, with one time sequence per column
}
#b) numberival application of the function we made in a)
m=5;K=50
sigma=0.3
a=c(0.4,0.5)
X=ex1317(m,K,sigma,a)$X
#plot.
matplot(-1:K,X,"l",xlab="Time",ylab="")
lines(0:K,0:K*0,lwd=3)
title("Second order autoregression")
#c) redo b) but with
a = c(0.4,-0.5)
#what differences do we see?
#what would happen if we instead had initial condition = 1?
#### 13.18
#a) Gonna run a single simulation using function from 13.17 a) up to time K.
# Then gonna look at the autocorrelation function of X. Parameters are the same as before.
m=1;K=10000
sigma=0.3
a=c(0.4,0.5)
X=ex1317(m,K,sigma,a)$X
acf(X)
#b) redo a) but with
a = c(0.4,-0.5)
# how can we link these plots to the ones we had in the exercise 13.17?
#### 13.19
#a) a general function which generate m sequences of a Gaussian ARMA(1,1) up to time K
#m = number of realizations
#a = AR-coefficient
#b = MA-coefficient
#sigma = standard deviation of the errors
#x_0 = initial value
ex1319=function(m,K,x_0,sigma,a,b){
X=matrix(x_0,K+1,m)
eps=matrix(rnorm(m*(K+1)),K+1,m)
for (k in 1:K){
se=sigma*(eps[k+1,]+b*eps[k,])
X[k+1,]=a*X[k,]+se
}
list(X=X)
#X = (K+1)xm matrix of simulations with one time sequence per column
}
#b) numerival application of function made in a) for given parameter values.
m=20;K=50
x_0=0
sigma=1
a=0.7;b=1
X=ex1319(m,K,x_0,sigma,a,b)$X
#Plot
matplot(0:K,X,"l",xlab="Time",ylab="")
lines(0:K,0:K*0+x_0,lwd=3)
title("First order autoregression and moving average")
#c) Single simulation using function from a) ut to time K, and look at auto-correlation
#function of X. Parameters are as in the earlier exercise
m=1;K=10000
x_0=0
sigma=1
a=0.7;b=1
X=ex1319(m,K,x_0,sigma,a,b)$X
acf(X)
#d) how does the plot differ from AR(1) acf?
#### 13.20
#a) a general function computing an independent sequence of K correlated pairs of
# standerd normals errors.
#rho = the correlation
ex1320=function(K,rho){
e1=rnorm(K)
e2=rho*e1+sqrt(1-rho**2)*rnorm(K)
eps=cbind(e1,e2)
list(eps=eps)
#epx= Kx2-matrix containing the correlated pairs
}
#b) use fundtion made in a) to make a function which generate a sequence og a Gaussian
# AR(1) pairs up to time K
#sigma = vector of standard deviations of the error processes
#rho = vector of error correlations of the error processes
#A = matrix of auto-regrissive coefficients
#Initial condition = 0
ex1320_b=function(K,sigma,A,rho){
eps=ex1320(K,rho)$eps
# intitial condition = 0
X=matrix(0,K+1,2)
#recursion.
for (k in 1:K) {
X[k+1,]=A%*%X[k,]+sigma*eps[k,]
}
list(X=X)
#(K+1)x2 matrix containing simulations
}
#c) numerical application of the a) and b) with two processes plotted jointly against time.
K=100;sigma=c(1,1)
A=matrix(c(0.8,0.5,0.0,0.8),2,2)
rho=0.5
X=ex1320_b(K,sigma,A,rho)$X
# Plot
matplot(0:K,X,"l",xlab="Time",ylab="")
title("First order bivariate auto-regression")
#which line belongs to which process?
#d) redo c) but with
A=matrix(c(0.8,0.5,-0.4,0.8),2,2)
#changes?
#### 13.21
#a) a general function which computes m sequences of log-normal rates of interest up
# to time K.
# We use AR(2) as a linear driver process
# xi = mean reate of interest
# sigma = standard deviation of the errors in the driver process
# a = vector of AR-coefficients.
# Initial values are the median of the distribution
ex1321=function(m,K,xi,sigma,a){
X=ex1317(m,K,sigma,a)$X
gamma=sqrt((1-a[2])/(1-a[1]**2-a[2]**2-a[2]-a[1]**2*a[2]+a[2]**3))
#standard deviation of driver process
sigma_x=gamma*sigma
#transformation:
r=xi*exp(-sigma_x**2/2+X)
list(r=r,sigma_x=sigma_x)
# r = (K+2)xm matrix of simulations. One sequence per column
# sigma_x = standard deviation of Xk when initial value is forgotten
}
#b) numerical application of function made in a)
m=5;K=50
xi=0.04;sigma=0.3
a=c(0.4,-0.5)
r=ex1321(m,K,xi,sigma,a)$r
#Plot
matplot(-1:K,r,"l",xlab="Time",ylab="rate of interest")
lines(-1:K,-1:K*0+r[1,1],lwd=3)
title("Second order auto-regression with log-normal transformation")
#printing the mean.
rmean=round(mean(r),digits=3)
print(rmean)
## 25.04.18 ## 13.22-13.23 + 13.26-13.28##
#13.22 is done by hand
#### 13.23
#a) Make a general dunction which generates a single sequence of interest rates up to time K.
# We do this under the Black-Karisinsky model we know from 13.22 starting in r_0.
# xi = mean rate of interest
# sigma = standard deviation
# a = auto-regression coefficient
# sigma_x = long-term standard deviation of driver process
# X = driver process (vector)
# rln = rates of interest (vector)
ex1323=function(K,r_0,xi,sigma,a){
sigma_x=sigma/sqrt(1-a**2)
x_0=log(r_0/xi)+sigma_x**2/2
X=rep(x_0,K+1)
for (k in 1:K) {
X[k+1]=a*X[k]+sigma*rnorm(1)
}
##Transformation:
#log-normal
r=xi*exp(-sigma_x**2/2+X)
#b) generates interest rates for same period under B-K and Gamma rates of interst.
#both computed from the same linear driver process
#gamma:
alpha=1/(exp(sigma_x**2)-1)
U=pnorm(X/sigma_x)
rgam=xi*qgamma(U,alpha)/alpha
#rln = log-normal rates of interests
#rgamma = gamma rates of interests
list(sigma_x=sigma_x,X=X,rln=r,rgam=rgam)
}
#c) numerical application of function made in a) and b) with plots
K=100;r_0=0.02
xi=0.04;sigma=0.3
a=0.7
o=ex1323(K,r_0,xi,sigma,a)
# Then the joint plot of the rates.
r=cbind(o$rln,o$rgam)
matplot(0:K,r,"l",xlab="Time",ylab="Rates of interest")
#to see which line belongs to the log-normal rates:
#lines(0:K,r[,1],col="blue")
title("Comparing log-normal and Gamma modelling of interest rates")
###13.26 and 13.27 done by hand
####13.28
# a) Make a general function which generates m simulated sequences up to time K
# of Wilkie inflation rates
# xi_i = initial rate
# sigma_i and a_i = Wilkie parameters
ex1328=function(m,K,xi_i,sigma_i,a_i){
X_i=matrix(0,K+1,m)
for (k in 1:K) {
X_i[k+1,]=a_i*X_i[k,]+sigma_i*rnorm(m)
}
I=(1+xi_i)*exp(X_i)-1
list(I=I)
#(K+1)xm matrix containging simulated rates of inflation, one sequence per column.
}
#b) numerical application of function made in a) with plots
m=10;K=50
xi_i=0.048;sigma_i=0.04
a_i=0.58
I=ex1328(m,K,xi_i,sigma_i,a_i)$I
#Plot:
matplot(0:K,I,'l',xlab="Years ahead",ylab="Rate of inflation")
lines(0:K,0:K*0,lwd=3)
title("The Wilkie rate of inflation")
#c) a general function using inflation from function in a) to simulate price levels
#parameters are as in the earlier exercises for this week
ex1328_c=function(m,K,xi_i,sigma_i,a_i){
I=ex1328(m,K,xi_i,sigma_i,a_i)$I
# Now the price level.
if(m>1){
Q=apply(1+I[1:K+1,],2,cumprod)
} else {
Q=cumprod(1+I[1:K+1,])
}
Q=rbind(rep(1,m),Q)
list(Q=Q)
#(K+1)xm matrix containing simulated price levels, one sequence per column
}
#d) numerical application of function made in c) with plots
m=100;K=50
xi_i=0.048;sigma_i=0.04
a_i=0.58
Q=ex1328_c(m,K,xi_i,sigma_i,a_i)$Q
#plot:
matplot(0:K,Q,'l',xlab=" Years ahead",ylab="Price level")
lines(0:K,0:K*0+Q[1,1])
title("The Wilkie price level")