UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical solution of a control problem in economics Han, Bing 1996

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1996-0338.pdf [ 3.61MB ]
Metadata
JSON: 831-1.0079781.json
JSON-LD: 831-1.0079781-ld.json
RDF/XML (Pretty): 831-1.0079781-rdf.xml
RDF/JSON: 831-1.0079781-rdf.json
Turtle: 831-1.0079781-turtle.txt
N-Triples: 831-1.0079781-rdf-ntriples.txt
Original Record: 831-1.0079781-source.json
Full Text
831-1.0079781-fulltext.txt
Citation
831-1.0079781.ris

Full Text

NUMERICAL SOLUTION OF A CONTROL PROBLEM IN ECONOMICS By Bing Han B . Sc. (Applied Mathematics/Industrial Management Engineering) Shanghai Jiao Tong University, 1991 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E DEGREE OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA June 1996 © Bing Han, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mci-the m a f ' The University of British Columbia Vancouver, Canada Date DE-6 (2/88) Abstract This thesis applies a numerical technique to a central bank control problem. The bank tries to control the inflation by influencing the nominal interest rates using the bank rates. The mathematical model is a two-dimensional singular control problem. The numerical method used is the Markov chain approximation method. We compute the optimal control and have some heuristic discussion on the convergence of the numerical method. Then we estimate the parameters in our model by the method of moments. We analyze the results from computing the control in the model with the estimated parameters and an economic interpretation is given. i i Table of Contents Abstract ii List of Tables v List of Figures vi Acknowledgement viii 1 Introduction 1 2 The Continuous Model and the Hamilton-Jacobi-Bellman Equation 4 3 Markov Chain Approximation 7 3.1 The numerical boundary and the numerical grid 7 3.2 The transition functions: case 1 10 3.3 The transition functions: case 2 13 3.4 The transition functions: case 3 14 3.5 The control problem for the chain 15 4 Dynamic Programming and Computation 17 4.1 The dynamic programming equation for the chain 17 4.2 Comparison of the dynamic programming equation and HJB equation 19 4.3 Preparation 19 4.4 Computational methods for dynamic programming equation 25 4.5 Discussion of convergence. . 27 iii 5 Test Results 32 5.1 Test 1: uncontrolled problem 33 5.2 Test2: diagonal drift matrix 34 5.3 Test3: different grid size 38 5.4 Test 4: heuristic discussion of the relationship between the deterministic differ-ential equations and the control 38 6 Data Fitting 44 6.1 Description of the data 44 6.2 Estimation method 45 6.3 Estimation Results 48 6.4 Control results and interpretation 50 7 Conclusion 56 Appendices 58 A Derivation of the Transition Probabilities for the Diffusion Case 58 B Verification of Consistency Conditions 60 C Cost Function of the Uncontrolled Problem: the Exact Solution 63 Bibliography 67 iv List o f Tables 5.1 Parameters for test 4 40 6.1 Four cases 47 6.2 Estimation results 49 v List of Figures 3.1 Numerical boundary 8 3.2 h-grid Gh and G+ 9 3.3 Grid points 11 3.4 Diffusion step 12 3.5 Control step 14 3.6 Reflection step 15 4.1 Elimination of the reflection states 20 4.2 Ordering of the state space 23 4.3 Rh 24 5.1 The region [-0.2,0.3] X [-0.2,0.3] 34 5.2 Cost function on the region [-0.2,0.3] X [-0.2,0.3] 35 5.3 Comparison of different grid sizes on horizontal lines 36 5.4 Comparison of different grid sizes on vertical lines 37 5.5 Control for uncorrelated problem 38 5.6 Comparison of different grid sizes 39 5.7 Comparison of free boundaries for different grid sizes 40 5.8 Control for test 4 41 5.9 Phase plane 42 6.1 Canadian treasury bill rates 44 6.2 Consumer price index 45 6.3 The Bank of Canada rates 45 6.4 The inflation rates 46 vi 6.5 Comparison of the four cases (/> = 0.1) 51 6.6 Comparison of different p's 52 6.7 Comparison of different /x's 53 6.8 Control on a bigger region 54 6.9 Comparison of different a's 55 vii Acknowledgement I am greatly indebted to my supervisor, Ulrich Haussmann, who has guided me through two and a half years of graduate studies in mathematics with patience. I am grateful to John Walsh and Brian Wetton for giving me suggestions on my thesis. I would also like to thank all the graduate students and the staff in the Mathematics Department for kindly helping me in many ways. viii Chapter 1 Introduction In recent years, using numerical methods to solve stochastic singular control problems has become an active area of research. Loosely speaking, we can think of a singular control as a sequence of small impulses, but the control cannot be represented as an integral of a "control rate", because the "rate" might not be a bounded function. The interest in this class of stochastic singular control problems has been increasing rapidly due to the fact that this type of problems has been frequently used in modeling the real world systems, especially economic systems (such as the one in this thesis). However, only limited analytical results are available for singular control problems. Even if the existence of an optimal control can be shown under some assumptions (such as in Haussmann and Suo [9] and [10]), usually we are not able to find that optimal control analytically. So people have started to look for a numerical approach. Generally, there are two methods of discretization. The first one is the discretization of a variational inequality that the value function satisfies. (In the language of control theory, the variational inequality is called Hamilton-Jacobi-Bellman (HJB) equation.) This method is an extension of the numerical schemes for deterministic differential equations. Examples in this category can be found in Tourin and Zariphopoulou [21], Bancora-Imbert, et al. [2] and Akian, et al. [1]. The second method is a direct discretization of the continuous stochastic process. Usually, we call this the Markov chain approximation. The approximate Markov chain can be constructed by purely probablistic methods which do not require the use of the analytical properties of the HJB equation. (Although sometimes the HJB equations can be suggestive of good approximate Markov chains.) This is one of the advantages of the second method over the first one since for many classes of problems, little is known about the analytical properties 1 Chapter 1. Introduction 2 of the H J B equations. In this category, there is an example by Fitzpatrick and Fleming [6], although they dealt with a regular control problem. A systematic and detailed explanation can be found in Kushner and Martins [16], and Kushner and Dupuis [15]. In this thesis, we apply the methods in [16] and [15] to a central bank control problem. In Chapter 2, we introduce the mathematical model for the central bank control problem. According to Binhammer [3] (cf. pp.309-310, pp.317-327), the Bank of Canada uses bank rates as one of its control factors to affect the short term interest rates so that the money supply growth is held in a certain range and thus the inflation is kept at some target level. We try to describe this behavior in terms of stochastic processes. Our model, based mainly on the work of Chiarolla and Haussmann [5], is a two-dimensional infinite horizon stochastic singular control problem. The stochastic processes of the nominal interest rate and the inflation rate are jointly dependent. The bank rate, which plays the role of control, is assumed to affect the interest rate directly and additively. The value function reflects the bank's desire to keep the interest rates and the inflation low over the long run and its reluctance to make large changes in the bank rate. Chapters 3 and 4 explain how the Markov chain approximation method in [16] and [15] is applied to the model introduced in Chapter 2. Chapter 3 deals with the construction of a Markov chain that approximates the continuous system. For the purpose of numerical computations, in Section 3.1, we modify our model a bit by restricting the process to a finite domain. Sections 3.2 through 3.4 obtain the transition probabilities for the approximate chain in different cases. The control problem for the approximate chain is established in Section 3.5. Methods for solving the discretized control problem are given in Chapter 4. In Section 4.1, the dynamic programming equation for the chain is derived. Section 4.2 compares the dynamic programming equation with the H J B equation. Sections 4.3 and 4.4 explain the iteration methods for computing the dynamic programming equation by iteration. Section 4.5 makes some observations on the convergence of the iteration methods and the convergence of the Markov chain approximation. Chapter 1. Introduction 3 The numerical results are included in Chapters 5 and 6. In Chapter 5, we test our ap-proximation method on the model using artificial parameters. There are three test problems and the results confirm the convergence comments made in Section 4.5. Chapter 6 analyzes our model using the data from Statistics Canada. Section 6.1 describes the data we use. Section 6.2 explains the method of moments that we use and the results are given in Section 6.3. Section 6.4 shows the control results by using the parameters derived in Section 6.3, and the economic interpretation is given. We conclude in Chapter 7. Chapter 2 The Continuous Model and the Hamilton-Jacobi-Bellman Equation The continuous model introduced in this chapter is from the work of Chiarolla and Haussmann [5]. Let (il, T, V) be a fixed complete probability space and (!Ft)t>o a given filtration. We suppose that the stochastic differential equations of nominal interest rate and inflation rate are dxx(t) = (ai + bnx1(t) + b12X2(t))dt + a1dW(t) + cdk(t) dx2(t) = (a2 + b2ix1(t) + b22x2(t))dt + a2dW(t) with initial values (2.1) (2.2) x i ( 0 - ) = xx x2(0~) = x2 where x\(i) and x2(t) represent the nominal interest rate and inflation rate at time t respectively. k(t) represents the bank rate at time t and the constant c measures the rate of influence that the change of bank rate will have on interest rate. W(t) = (wi(t) w2(t))' is an 7£ 2-vafued standard Brownian motion with w\(t) and w2(t) independent of each other, representing uncertainty in the economy. a\ = ( C T U O\2) and a2 = (a2\ a22) are volatility coefficients. 'a\, 6 n , a2, b2\, b22 are constants. Let / X(t) xi(t) \ * 2 ( « ) x2 i \ o 2 ] ( , K(t) = ck(t) 0 B(X(t)) 1 a-i + buXi(t) + b12x2(t) ^  K a 2 + b21xi(t) + b22x2(t) / 4 Chapter 2. The Continuous Model and the Hamilton-Jacobi-Bellman Equation 5 then (2.1) and (2.2) can be rewritten in a short form dX(t) = B(X(t))dt + adW(t) + dK{t) (2.3) X(0~) = X. (2.4) The set U(X) of admissible controls is the collection of adapted stochastic processes k{t) which are right continuous with limits on the left, having locally bounded variation and k(0~) = 0. The process X(t) is thus right continuous with left hand limit and (2.3) and (2.4) are equivalent to: X(t) = X+ [* B(X(s))ds + o ftdW(s) + K(t) Jo Jo for t > 0. The cost function corresponding to the control k is defined as UO O 1 POO -j -(p.x21(t) + x22(t))e-ptdt + aj e-^dvk(t)j (2.5) where Ex denotes expectation given the initial condition X = (xi,x2) and Vk{t) = var(& : [ 0 ~ , 2 ] ) . a and p are weights with a > 0 and fi > 0. p is a given positive discount factor. The problem is to minimize J over all admissible controls and to find the value function V(X)= inf Jx(k). (2.6) U(X) We are here facing an infinite horizon singular control problem which is similar to the one in Chiarolla and Haussmann [5], although they used a finite horizon. Compared to their optimal cost function, ours is time independent. That is why in (2.6) V is only a function of X. From their results, we expect that our V(X) has weak second derivatives in £°°(7£2) and satisfies the following variational inequality min^CV - pV+^(p,xl + xl),mm{a-VXl,a + VXl}^ = 0 (2.7) where 1 2 d2V CV = (ai + buxx + b12x2)VXl + (a 2 + b2ix1 + b22x2)VX2 + - aii Q X . Q X (2-8) i,j-l 1 3 Chapter 2. The Continuous Model and the Hamilton-Jacobi-Bellman Equation 6 and 2 k=l so if a — (a{j), then a - aa . (2.9) In the following chapters we are going to approach this singular control problem from a numerical perspective. Chapter 3 M a r k o v C h a i n Approx imat ion In this chapter we are going to find a discrete Markov chain that approximates the system (2.1) and (2.2) and the associated discretized control problem for which the desired computation can be carried out. The approximating Markov chain is chosen such that certain "local" properties of the approximating chain are "similar" to those of the original controlled process. We follow Kushner and Martins [16] and Kushner and Dupuis [15]. 3.1 T h e numerical boundary and the numerical grid. The state space for the process (2.1) and (2.2) is on the whole plane. However, for the purpose of numerical computation, we need to constrain the process to take values in a finite domain which should be big enough to retain the essential features of the process and the corresponding control problem. The boundary of the finite domain thus defined is called the numerical boundary because it is not a real boundary. The finite domain we use is a two dimensional box G = [ A n , A21] * [^ .12,^ .22] as shown in Figure 3.1. (Nb. Restricting the domain to a suitably large rectangle has little effect on the result, cf. Figure 6.8.) Wi th the numerical boundary, we must add a boundary condition. A direct and intuitive choice is to make the boundary reflecting. So (2.1) is revised to the following: cfai(t) = (a1 + bux1(t) + b12x2(t))dt + o-1dW(t) + cdk(t) + dY1(t)-dU1(t) (3.1) dx2(t) = (a2 + b21x1(t) + b22x2(t)dt + a2dW(t) + dY2(t)-dU2(t) where l i ( - ) (i = 1,2) and £/;(•) (i = 1,2) are reflection terms which are nondecreasing. !;(•) (i = 1,2) can increase only when Xi(t) = An (i = 1,2) and Ui(-) (i = 1,2) can increase only when X{(t) = A2i(i = 1,2). This is illustrated in Figure 3.1. The other symbols in (3.1) are the 7 Chapter 3. Markov Chain Approximation 8 x 2 ( A 11 , A 22 ) u2 ( A 21 , A 22) G U l ( A U i A 1 2 ) Y 2 ( A 2 1 > A 1 2 ) Figure 3.1: Numerical boundary same as in (2.1). Here we have assumed that the reflection direction is perpendicular to the boundary. Note that at the corner of the rectangular boundary, the reflection direction can be any combination of !;(•) (i = 1,2) and [/,•(•) (i = 1,2) as long as it points into G. Since the discrete Markov chain will be defined on a discrete state space, first we need to discretize the region G. This is shown in Figure 3.2. Let h be an approximation parameter {h is a scalar here, but it could be vector valued, that is, the grid spacing could depend on the direction, cf. Kushner and Dupuis [15], Section 5.3.1), then define the /i-grid Gh on G by Gh = { ( A n + m/i , Ai2 + lh) : 0 < m < (A2\ - Au)/h, 0 < I < (A22-Ai2)/h} where we assume, without loss of generality, that (A21 — An) and (A22 — A\2) are integral multiples of h. Define the "extended /i-grid" G% = {(An + mh,A12 + Ih) : - 1 < m < (A21 - Alx)/h + 1, - 1 < / < (A22 — A\2)/h + 1}. Define dGh = Gh fl dG, the grid points on the boundary of G and the "reflecting boundary" dG~fc = G^ — Gh- The extended grid G^ and dG~fc are introduced to simpUfy the exposition. Chapter 3. Markov Chain Approximation 9 (Au,A22) A A (An,A12) x2 A A A A A A A A A o o o o o o o o o o o o o o o o o o o o o o o o o o G. o o o o o o o o o o o o o o o o o o o o o o o o (A2i,A22) (A21, A i 2 ) <-dGh dGi Figure 3.2: h-grid and G£ Wi th all the preparation above, we are ready to construct a Markov chain that "behaves like" (3.1), a controlled diffusion process with reflection on the boundary. Intuitively, the control term k(-) and the reflection terms Yi(-)(i = 1, 2) and Ui(-)(i = 1,2) act either "instantaneously" or "impulsively". So we can "split" the process into three cases: 1. Diffusion without control inside the domain. 2. Pure control inside the domain. 3. Reflection at the boundary dG%. The chain and the control will be chosen so that only one of the three cases happens at each step. In the following sections, £^ (n < 0 0 ) is the discrete approximate Markov chain on G\. Chapter 3. Markov Chain Approximation 10 3.2 The transition functions: case 1. This is the case when = X € Gh and no control is applied. Following the method in Kushner and Dupuis [15], the transition functions are chosen such that the first two moments of A££ = - ^ are close to those of / Bdt + f adW over a small time interval. To be precise, the following "local consistency" conditions have to be satisfied: EXtnA& = B{X)Ath{X) + o(Ath(X)), (3.2) Ex<n[A^n - £*iBA£*][Atf - EXinAtt}' = aAth(X) + o(Ath(X)) where a is the same as defined in (2.9) and E\ n denotes the conditional expectation given {g,i<n,£ = X}. Following the examples in Kushner and Dupuis [15], we constructed the transition functions by using a finite difference approximation of the derivatives in the uncontrolled problem's Hamilton-Jacobi-Bellman(HJB) equation: _ CV + pV - h^nxl + x\) = 0 (3.4) where CV is defined by (2.8). Using a one sided scheme to approximate the first derivatives gives, . dv(x) dxi \y(x + eih)-v(x)]/h UBi(x)>o (i = l ,2 ) (3.5) [(V(X) - V(X - eih)] /h if Bi(X) < 0 where (Bi(X),B2(X)) = B(X)(as in (2.4)) and e;(i = 1,2) is the i t h standard basis element in n2. The difference scheme approximating the second derivatives is the following. If i = j(i,j = 1,2), d2V(X) V(X + eih) + V(X - eih) - 2V(X) dxidxi h2 This is a standard central difference scheme for second derivatives. (3.6) Chapter 3. Markov Chain Approximation 11 If i ^ j and > 0 = 1,2), d2V(X) 2V(X) + V(X + ejh + ejh) + V(X - e{h - e3h) dx 'Ox' V{X + ejh) + V(X - e^ fe) + F ( X + ejh) + F ( X - ejh) 2h2 To make the above approximation clearer, the right hand part can be rewritten as 1 f F ( X + ejh + e3h) - V(X + ejh) _ V(X + gh) - V(X) 2 I h h •}/h (3.7) F ( X ) - V(X - e,-fe) _ V(X - ejh) - V(X - ejh - ejh)' h h So the idea is to approximate | ^ first, then approximate the mixed second derivative. The grid points involved in the above scheme are shown in Figure 3.3(a). Note that the diagonal which has points on it is in the "forward direction". Hence this scheme for mixed second derivatives shares the same principle as the one sided first derivative scheme. If i ^ j and a tj < 0 ( i , j = 1,2), similarly, we have d2V(X) 2V(X) + V(X + e{h - ejh) + V(X - ejh + ejh) dxidxj 2h2 so o\ V(X + ejh) + V(X - ejh) + V(X + ejh) + V(X - ejh) 1 ' ; + 2h2 The grid points used is shown in Figure 3.3(b). Note that the diagonal which has points on it is in the "backward" direction. (a) (b) Figure 3.3: Grid points We assume that at each step the approximate chain only moves to the adjacent grid points as shown in Figure 3.4 (so there are eight possible positions each time). After some calculation Chapter 3. Markov Chain Approximation 12 (cf. Appendix A) and recalling that a12 = a21, we get the following transition probabilities. Notation ph((xi, x2), (j/ 1 ? y2)) stands for the transition probability from point (21,2:2) to point (2/1,2/2) and X = (x1,x2). time t+ A t (X) h— h —H Figure 3.4: Diffusion step ph((x1,x2),(x1 + h,x2)) = ph((xi,x2),(x1 - h,x2)) = ph((x1,x2),(x1,x2 + h)) ph((x1,x2),(x1,x2 - h)) + hB+(X) af-\^ + hB;(X) IQ\X) IQ\X) af-^ + hBt(X)\/Q\X) 022 |«2l| + hB2(X) IQ\X) Ph({Xl,x2), (xt + h,x2 + h)) = \a+2/Qh(X) ph((xux2),(Xl - h , x 2 - h)) = -a+2/Qh(X) ph((xux2),(Xl -h,x2 + h)) = ^2/Qh(X) (3.9) (3.10) (3.11) (3.12) (3.13) (3.14) (3.15) Chapter 3. Markov Chain Approximation 13 ph((x1,x2),(x1 + h,x2- h)) = \a^2IQh(X) Qh(X) = a n + a 2 2 - ± ( | a 2 1 | + |a 1 2 | ) + fc(l*iPOI + I ^ W I ) . (3.16) where the normahzation factor (to make the ph's sum to 1) is (3.17) To make the probabilities well defined, that is, to make the probabilities nonnegative and the denominators in the above formulas nonzero, we need a n > K 2 I «22 > |«2i| (3.18) If this restriction does not hold, we must look for other consistent transition probabilities (cf. Kushner and Dupuis [15], pp.109-110). Fortunately, the parameters we get for the model satisfy this condition (cf. Chapter 6). If we define the interpolation interval Ath(X) h2 (3.19) Qh(X) then the chain defined by (3.9)-(3.16) is locally consistent with the uncontrolled diffusion process in a sense described by (3.2) and (3.3). Actually, it can be shown (cf. Appendix B) that Ehx„A£ = B{X)At\X) (3.20) E X,n aAth(X) +o(Ath(X)) (3.21) 3.3 The transition functions: case 2. This is the case when £^ = X 6 Gh and control is exerted. Since the control is "impulsive", we let Ath(X) = 0. Referring to (3.1), we can see that the control is one dimensional, so there are only two possible directions of the control: either increasing or decreasing along xi-axis. We let Ak% denote the control action along x\ direction at step n. For programming simplicity, we restrict Akl\ so that the process only moves to the adjacent grid points. See Figure 3.5. If Chapter 3. Markov Chain Approximation 14 control is exerted at X, then the process moves to the grid point either to the right or to the left of X. The restriction on the size of the control at each step that we have used here will yield the same limit results as using arbitrary values for Ak[\. (cf. Kushner and Martins [16], Section 5.) time t Figure 3.5: Control step We can express the above idea in terms of conditional probabilities. That is ph((x1,x2),(x1 - h,x2)\Akl) = -h) = l ph((x1,x2),(x1 + h,x2)\Ak^ = h) = l. We will see later that the control defined here give rise to a natural numerical procedure. In fact, it is very similar to a finite difference approximation of (2.7). 3.4 The transition functions: case 3. This is the case when the approximate chain reaches the artificial boundary dG^. The numerical step we are going to use is a simple and straightforward approximation of the reflection terms in (3.1). It takes X = ^ €. dG^ to the nearest point on Gh as shown in the figure 3.6. If we let A Y ^ (i = 1,2) and AU^^i = 1,2) denote the reflection at step n, then we have Chapter 3. Markov Chain Approximation AY£{ = h(i = 1, 2) and Af/£- = /i(i =1,2) when necessary. 15 Figure 3.6: Reflection step Let AY// = (AY//A, AY£2) and AU^ = ( A [ / £ 1 ; A f / ^ 2 ) and combine section 3.3 and section 3.4, then we can define Ynh = £ AY/1 Vhn = J2 AU? khn = ]T Ak/ i=0 i=0 i=0 3.5 The control problem for the chain. We say that a control policy kh is admissible if it preserves the Markov property P{tn+i = Y\g,Ak$,i< = X} = ph(X,Y\Akhn). Let At% = Ath(£%) if n is a pure diffusion step, and set At1^ = 0 otherwise. Define i£ = Y^i^o1 Atf. Note that {Afc^;} is a discrete time process here and is known at time i£. Wi th initial condition $ = X, the cost function for the approximating chain under control policy kh is { oo i oo "j £ 5 e~"'"(M(£n,i) 2 + ( t f , 2 ) 2 ) A i * ( t f ) + « £ ^ p i " | A ^ | (3.22) ra=0 n=0 J where (i = 1,2) represents the component of ^ and the objective is to minimize j \ Chapter 3. Markov Chain Approximation 16 over all the admissible kh,s and to find the value function Vh(X) = mfJx(kh). (3.23) (3.22) and (3.23) are discretized versions of (2.5) and (2.6). Next we are going to discuss methods for dealing with this discrete optimization problem. Chapter 4 Dynamic Programming and Computation In this chapter, a dynamic programming equation for (3.22) and (3.23) will be derived and numerical methods for solving the equation will be suggested. Also there is discussion on convergence. 4.1 The dynamic programming equation for the chain. Assume for now that X S Gh, and we will deal with the case when X is on the boundary (that is, when X 6 dG^) later on. We now proceed to obtain a functional equation for VH(X) in (3.23). At current state X, we need to decide whether to apply the control or not. If we do not apply the control, whether or not the decision is optimal, then there is an immediately realized cost ±(nxl + x22)Ath(X) (4.1) as well as costs to be realized in the future. In order to get the optimal value function (least cost), we will have to use an optimal policy from the next step on, from whatever position the system finds itself in at that time. This is called "the principle of optimality". In this spirit, the mean value of the future cost, discounted to the present, is e-PATH^E(Vh(Cn+i)\U = X) = e - P M H ^ E X V H { ^ ) = e"'*** W £ Vh{Y)p\X, Y) (4.2) Y where the sum is over all the states that the chain can reach in one step from X. By (4.1) and (4.2), the total cost for not applying any control is c - , A t * ( x ) J2PH(X, Y)V\Y) + \(pix\ + xl)At\X) (4.3) Y L where X — (xi,x2). 17 Chapter 4. Dynamic Programming and Computation 18 If control is exerted, we have two choices, either move the chain to the left or to the right. If we move it to the right, then the total cost is the cost of the control plus optimal cost of the right point, that is ah + Vh(Xl + h,x2) (4.4) The optimal cost function Vh(x\ + h, x2) is not discounted to the present as in the diffusion case, because the control is "impulsive", which means no time has elapsed during this step. Similar argument will lead to the cost of moving to the left, so it is ah + Vh(x1-h,x2) (4.5) The control direction we choose should be the one that has least cost, so the cost for exerting control is min [Vh(xi + h, x2) + ah, Vh(x1 - h, x2) + ah} (4.6) The final optimal decision to be taken is the decision which attains the minimum of the costs over the two possibilities, which are: (1) exerting control, (2) no control and then acting optimally in the future. Hence, the optimal cost must satisfy the equation V\X) = min { e - ^ W ZY ph(X, Y)Vh(Y) + \(iix\ + x22)Ath{X), min [vh(xi + h, x2) + ah, Vh{xx - h, x2) + ah } = min | e - ^ P 0 £ y p\X, Y)Vh{Y) + \(px\ + x2)Ath(X), Vh(x1 + h, x2) + ah, Vh(x1 - h, x2) + ah} . This is called the dynamic programming equation and it is true for all X £ Gh-When X is on the numerical boundary, that is, when X £ dG^, since the reflection step costs nothing and requires no time, by "the principle of optimality", the cost function is Vh{X) = ExVh(£1) = YJph{X,Y)Vh{Y) (4.8) Y where the sum is over all the points Y that the chain is reflected to from X. In our case, the reflection step takes the chain to the nearest grid point on Gh, so there is only one Y in (4.8) (4.7) Chapter 4. Dynamic Programming and Computation 19 and (4.8) becomes Vh(X) = Vh(the nearest point to X on Gh) (4.9) This says that when X is on the numerical boundary, we take the value function Vh(X) the same as the value function of the nearest grid point. 4.2 Comparison of the dynamic programming equation and HJB equation. By moving Vh(X) = Vh(x\, x2) to the other side in equation (4.7), it becomes min \e-^ihW £ y ph(X, Y)Vh(Y) + \{u.x\ + x22)Ath{X) - Vh(X), min Vh(x1 + h, x2) - Vh(x1, x2) + ah, Vh(x1 - h, x2) - Vh(x1,x2) + ah} } = 0 (4.10) The first part in the braces above can be considered as a "diffusion part". If it is divided by Ath(X), after a long calculation, it can be seen that it is a finite difference approximation to the first part in (2.7), which is the same as (3.4)-the uncontrolled problem's HJB equation. It is not surprising because we derived the transition functions for the diffusion step by discretizing (3.4). If the second part in the braces above is divided by h, then it is very straightforward to see that 'Vh(Xl + h,x2) - Vh(x!,x2) + ah Vh(Xl - h,x2) - Vh(xux2) + ah\ min h h J nmmia + V^.a-V*) which is very similar to the second part in the HJB equation (2.7) if the superscript h is dropped. Moreover, (4.9) implies that ^n~\dGh = 0 where n is the unit outward normal. 4.3 Preparation. Before we go into the details of the methods for solving (4.7) and (4.9), we need to make some changes and introduce some notation. Chapter 4. Dynamic Programming and Computation 20 (4.7) and (4.9) are equations of X £ G%. For programming simplicity, we substitute (4.9) into (4.7) where appropriate and thus ehminate the reflection states (that is, all X 6 9G^). We are going to illustrate the method by an example, cf. Figure 4.1, (a). XQ is a reflection (a) (b) Figure 4.1: EUmination of the reflection states state in the middle of the boundary. By (4.9), we have V\X0) = Vh(X2). (4.11) We want to substitute (4.11) into (4.7) where Vh(Xo) shows up in (4.7). According to the definition of the transition probabilities in Chapter 3, Xo has connection with X\, X2 and X 3 . So there are three equations in (4.7) that have Vh(X0) in them, that is, Vh(X\), Vh(X2) and Vh(X3). Check Vh(Xi) first. The chain can go from X\ to Xo and X2. By (4.7), we can write V\X{) = min {e->**l*i) E r ^ 0 , x 2 ) / ( ^ i ^ ) ^ ( ^ ) + 6 - ^ f e ( ^ ) ( / ( x 1 , X o ) y f e ( x 0 ) + / ( x 1 , x 2 ) y f e ( x 2 ) ) + £ ( M 4 i + 42)A*'l(Xi), Vh(xhl + h, xh2) + ah, Vh(xhl - h, xh2) + ah} with X i = (£1,1, xii2). By (4.11), the underhned part becomes (ph(XuX0) + ph(X1,X2))Vh(X2). Chapter 4. Dynamic Programming and Computation 21 So we can define p^(X1,X2) = ph(X1,X2)+ph(X1,XQ). .(4.13) Thus Vh(X0) is ehminated from the right hand side of Vh(X1) in (4.7). Similarly, we can define p~t(X3,X2)=ph(X3,X2) + ph(X3,X0) (4.14) so that Vh(Xo) is ehminated from (4.7). As to Vh(X2), the chain might move from X2 to Xo , but it does not stay at X2. By (4.7), we can write V\X2) = min { e - ^ t " ( X 2 ) ph{x^ Y)Vh(Y) +e-"Ath(x^ph(X2,X0)Vh(X0) t 1 (4-15) + l(»x22A + xl2)Ath(X2), Vh{x2tl + h,x2>2) + ah,Vh(X0) + ah} with X2 = ( x 2 , i , Z 2 , 2 ) . So Vh(X0) appears in both diffusion and control terms. By (4.11), the first underlined part in (4.15) becomes p\X2,X0)Vh(X2). So in a diffusion step, we define ?(X2,X2)=ph(X2,X0). (4.16) Similarly, we substitute the second underlined part in (4.15) by (4.11), and define in a control step p*(X2,X2\Ak* = -h) = l. (4.17) By (4.16) and (4.17), Vh(X0) is eliminated from the right hand side of Vh(X2) in (4.7). We have to be more careful when we define ph(X2,X2) in a diffusion step when X2 is a corner point as shown in Figure 4.1, (b). As we can see, XQ communicates with X'Q, XQ and Chapter 4. Dynamic Programming and Computation 22 Xg'. By (4.7), Vh{X2) = min { e - ' A t t ^ ) E F « ^ , l « i r } / ( ^ , F ) ^ ( F ) +e-*thV*){jfi{X2,X'Q)V\X'Q) + ph{X2,X>>)Vh{X>>) +ph(X2,Xg')Vh(X[l')) (4.18) + 1(^ 1,1 + a ;^ , 2 )Ai / l (X 2 ) , F h ( a ; 2 ) 1 + / i , z 2 , 2 ) + ah, Vh{X'0) + ah} . By (4.9), we have V\X'Q) = Vh(XH) = Vh(X'f) = Vh(X2). So the first underUned part in (4.18) can be written as (ph(X2, X'Q) + p\X2,X'Q) + p\X2,XZ'))Vh(X2). Hence in a diffusion step, we define ?(X2,X2) = ph(X2,X^) + ph(X2,Xg) + ph(X2,X[;'). (4.19) The definition for the control step is the same as before. Similarly, we can eliminate other points on dG^. By this elimination, we will actually solve (4.7) with the newly defined transition probabilities ph,s. (Nb. ph = ph when boundary points are not involved.) To calculate V(X) for each state X, we want an ordering of Gh- The ordering we use is shown in Figure 4.2. Assume that there are N intervals in the x\ direction and M intervals in the x2 direction. We start from the left bottom corner node (0) , go horizontally all the way to the rightmost node (N), then jump to the leftmost node (N + 1) on the next level up, all the way to the rightmost node (2(N + 1) — 1). Repeat this procedure until we reach the top right corner node ( ( M + 1)(N + 1) - 1). Wi th this ordering and the elimination of dG%, we can express (4.7) in vector form. Define a matrix Rh = {r\X,Y),X,Y € Gh) by rh(X,Y) = e~pAth(x)ph(X,Y). Because we have ( M + l)(iV + 1) states in dGh, the dimension of Rh is (M + 1)(N + 1) by (M + 1)(N + 1). The Chapter 4. Dynamic Programming and Computation 23 ( A N A 2 2 ) ( A 2 1 J A 2 2 ) node(M+l)(N+l)-l x l node 2(N+1)-1 ' ( A 2 l , A 1 2 ) node N Figure 4.2: Ordering of the state space we can see. form of Rh is shown in Figure 4.3. Each block represents one level of the grid. As Rh is a nine-diagonal matrix. Then we can write (4.7) as Vh = min {RhVh + C \ mm{MrVh + ah x I , MfVh + ah x 1}} (4.20) where the minimum is taken elementwise. Vh = {Vh(X),X £ Gh} is a 1 by ( M + 1)(N + 1) vector. And Ch = {\{nx\ + xl)Ath(X),(xl,x2) = X e Gh} is a 1 by ( M + l)(N + 1) vector as well. Mr and Mi are defined as the following: Mr U u Mi (M+i)blocks L (4.21) J (M+i)blocks Chapter 4. Dynamic Programming and Computation 24 where U 0 1 0 L = '•• 1 1 -J (JV+l)x(JV+l) 1 is a 1 by ( M + 1)(^V + 1) vector with 1 as its elements. 1 0 (4.22) (N+l)x(JV+l) Chapter 4. Dynamic Programming and Computation 25 4.4 Computational methods for dynamic programming equation. There are two types of classical methods: the approximation in policy space method and the approximation in value space method. We use both methods for our model. i) Approximation in value space. We tried two methods under this category: the Jacobi iteration and the Gauss-Seidel iter-ation. Algorithm for Jacobi Iteration: 1) Choose an initial value function V0H(X). 2) The (n + l)st iteration is Vnh+1(X) = min { e - p A t * ( X ) Zy ?(X, Y)VNH(Y) + \{px\ + z | ) A ^ ( X ) , min [ v ^ z i + h, x2) + ah, Vn(xx - h, x2) + ah] } . Note, if X is on the left boundary of Gh, that is, when x\ = An (cf. Figure 3.1), then in a control step, (2:1 — ft, x2) will be on dGf. By (4.9), we use V^(the nearest point to (x\ — h,x2) on Gh) instead. Similar argument can be applied when x\ — A2\. 3) If maxx | V r a / l + 1 ( X ) - V£(X)\ < c, stop and take V£+1(X) as the result. Otherwise repeat step 2. (e is the stopping criteria set by the user and the maximum is over all the states in the state space.) Algorithm for Gauss-Seidel Iteration 1) Choose an initial value function VQ(X). 2) The (n + l)st iteration is VNH+1(X) = min { e - ^ P D (ZY<X?(X ,Y)VNH+1(Y) + ZY>X?(X,Y)VZ(Y)) + \(lixl + xl)Ath(X),mm [v£(xi + h,x2) + aft, V ^ z a - h,x2) + ah] } (4.24) where " Y < X " means that point " Y " is ordered before " X " and " Y > X " has similar meaning. If X is on the left or right boundary of Gh, we need to treat it the same as in Jacobi method. Chapter 4. Dynamic Programming and Computation 26 3) If maxx \V£+1(X) - V£(X)\ < e, stop and take V£+1(X) as the result. Otherwise repeat step 2. (The meaning of e and maximum is the same as in Jabobi iteration.) ii) Approximation in Policy Space. In our model, there are only three choices for the control: no control, force to the left and force to the right. We use k = 0, k = — h and k — h respectively to denote the above three cases. Algorithm 1) Choose initial control for each grid point. For simplicity, We start from ko = 0 at each point, then solve V0\X) = e-»AthW J2?(X, Y)V0h(Y) + \(px\ + x22)Ath(X) Y 1 for V0h(X). 2) At the (n -f l)st iteration, compute the control first. Let k^+1 = 0, if e-^hWYy(X,Y)V£(Y) + \{iix\ + x2)Ath(X) Y < min ( V ^ ( z i + h, x2) + ah, V ^ x i - h, x2) + ah) . If the above inequality is not true, then let ki\+1 = h, when V£(xi + h, x2) + ah> V„(xi - h, x2) + ah and, let kl\+1 = —h otherwise. Wi th the control just defined, solve the linear system e - ^ ) E y JJ ( X ) y ) V * + i ( y ) + + ^ ( j r ) if = 0 V£+i(X) = I V*+1(Xl + h, x2) + ah if khn+1 = h (4.25) Vn+1(x-i - h,x2) + ah if kl\+1 = -h Note that similarly, when X is on the left or right boundary of G, we need to treat it specially as in Jacobi method. Chapter 4. Dynamic Programming and Computation 27 3) If maxx \V^+1(X) - V£(X)\ < e, stop and take V£+1(X) as the result. Otherwise repeat step 2. (The meaning of e and maximum is the same in Jacobi iteration.) The above algorithm requires the solution of a linear system to obtain V£+1(X) in step 2. Since the system is very big, the usual way to solve it is to use iterative methods, such as Jacobi or Gauss-Seidel iteration methods. 4.5 Discussion of convergence. There are two types of convergence we need to consider. The first one is the convergence of the iteration method that is used to find Vh(X) for each h. The second one is the convergence of Vh(X) to the exact value function V(X). As to the convergence of the iteration method, it depends very much on the iteration matrices. So first we are going to explain what is an iteration matrix. Let us start from the method of Jacobi iteration in value space. As in (4.20) We can write (4.23) in vector form, V*+1 = min {RhV£ + Ch, m i n { M r V * + ah x 1, Mtv£ + ah x 1}} ^ = min {RhV£ + Ch, MrV% + ah x 2, MtV* + ah x l } where the minimum is taken componentwise. To explain the concept of iteration matrix, we want to write (4.26) as K + i = Rh(kn)V£ + Ch{kn) (n = 1,2 ...) (4.27) where Rh(kn) and Ch(kn) are composed in the following way. Suppose we consider ( V ^ + 1 ) ; . The i t h row of Rh(kn) is chosen from the ith row of Rh (cf. Figure 4.3), M ; or Mr (cf. (4.21) and (4.22)) and the element of Ch(kn) is chosen from the element of Ch or the number ah, depending on where the following minimum is achieved: min{ i t h row of Rh • Vh + i t h row of C \ i t h row of Mr-Vh + ah, i t h row of Mi • Vh + ah] where '•' is a dot product of a row vector and a column vector. So Rh{kn) is a ( M + 1)(-/V + 1) Chapter 4. Dynamic Programming and Computation 28 by ( M + 1)(JV + 1) matrix (M and N are the same as in Figure 4.2), and we call this Rh(kn) the iteration matrix. For each n, Rh(kn) might be different. For each row of Rh(kn), as we have seen, there are three choices, and there are ( M + l)(iV + 1) rows of Rh{kn). So theoretically there are 3(^+1)^+1) choices for Rh(kn) (n = 1,2...). This Rh(kn) will also appear in the method of Gauss-Seidel iteration in value space and the method of iteration in policy space as shown below. For the method of Gauss-Seidel iteration in value space, we can write (4.24) in vector form as well. By an argument similar to the one used for the Jacobi method, we have V}+1 = L\kn)V^+1 + Uh(kn)v£ + Ch(kn) where Lh(kn) is the strict lower triangular matrix (not including the main diagonal) of Rh(kn), and Uh(kn) is the upper triangular matrix including the main diagonal. Similarly, for the method of iteration in policy space, we can write (4.25) in vector form vf+i = Rh(kn+1)v£+1 + c\kn+1). The usual convergence argument (such as the one in Kushner and Dupuis [15], pp.154-157) needs the contraction property of Rh(kn), that is (Rh(kn)y ^ 0 as i -»• 00. We only know that Rh(kn = 0) = Rh as defined in (4.20) has this contraction property. The reason is stated below. According to the definition of Rh, the sum of the absolute values of the elements in a row is £ \rh(X,Y)\ = £ e - ^ W ? ( X , Y ) = e - < ^ W £ ?(X,Y) YeGh YeGh YeGh where Y^YeGh Ph(X, Y) is the probability that X moves to all the state in the state space which gives 1. So we have • £ \rh(X,Y)\ = e-pAthW <1. YeGh Chapter 4. Dynamic Programming and Computation 29 This is true for every row, so the norm of R is less than 1, that is H ^ H o o = max £ \rh(X,Y)\ < 1. X Yeoh Let |A| denote the largest magnitude of the eigenvalue of Rh. Then |A| < H i f ^ H o o < 1, (cf. Noble and Daniel [18], p.268). Hence (Rhy -»• 0 as i -> oo (cf. Noble and Daniel [18], p.378). For arbitrary Rh(kn), we can see from the construction of Rh{kn) that some of the rows of Rh are replaced either by the corresponding rows of Mr or the corresponding rows of M ; . From (4.21) and (4.22), we see that each row in Mr or Mi has all entries zero except one which is 1, but that makes | | i2 / l (A; n ) | | 0 0 = 1. So the above argument for Rh does not apply any more. Also Rh(kn) is non-symmetric. Besides that, as we mentioned before, the matrix form of Rh(kn) has a huge amount of choice in the minimization. A l l these factors make the study of the properties of Rh(kn), and correspondingly the study of the convergence of the iteration methods, very hard. Fortunately, our numerical tests seem to converge. In fact, pathological non-convergence could only occur if the process does not spend an infinite amount of time with k = 0, i.e. away from the "free boundary". Our discussion of the convergence of Vh(X) to the exact value function V(X) is based on Kushner and Martins [16], Section 5, and Kushner and Dupuis [15], Section 11.2. We define kh(t)= £ Akhn (4.28) » : * n + l < * where is defined in Section 3.5 and Ak% is defined in Section 3.4. Now let us examine some properties of our model. • The process defined in (3.1) has continuous drift and volatility because the drift is linear and the volatility is constant. And we have a continuous utility function \{nx\ + x2). • sup/j Ex\kh(t)\ < oo for each t where | • | denotes total variation. The reason is the following. First, we show that sup^ J | ( 0 ) (0 stands for no control) as defined in (3.22) is Chapter 4. Dynamic Programming and Computation 30 finite. By (3.22), we have 4(0) = EX | £ \e-ptH^n,i)2 + ( t f ,2 ) 2 ) A **(tf )} • ( 4 - 2 9 ) Since the chain is defined in the finite box G, we have Md , i ) 2 + ( ^ , 2 ) 2 < C i where C\ is a positive constant which is independent of h. This gives 4(0) < \c.Ex jf>-**Ai*(tf)J . (4.30) As A i ' l ( ^ ) —y 0, the sum in (4.30) is an approximation to the integral poo 1= / e-"*dt. Jo The positivity of the p in our model guarantee that I < oo. Hence oo J2e-pt"Ath(^)<C2 n=0 where C2 is a positive constant independent of h. So 4(0) < \cxc2. Because C\ and C2 are independent of h, we have sup^ J^-(O) < oo. By (3.23), BViTpVh(X) < sup Jx(0) < oo. (4.31) h h As everything in (3.22) is positive, by (4.31), we can say that suphExi:^=oe-<\Akhn\<oo s u P / l Ex < t e" p t" IAfc£| < oo n- | - l — =>• e _ p i sup^ E x « < 0 0 because e_"*» > e~pt =* c-"* sup f c £x | f c f c ( i)l < because E n A x , |A*£| = |*h(*)|by (4.28) n + l — =>• sup^ Ex\kh(t)\ < oo because e~pt is a fixed number. So we have shown that sup^ Ex\kh(t)\ < oo. • i^n+i — '• n} I S a uniformly integrable family of random variables. This is because at each step our approximate control only takes value 0 or h. Chapter 4. Dynamic Programming and Computation 31 • The uncontrolled problem has a weak sense unique solution for each initial condition, (cf. [15],p.27.) • For each 6 > 0, there is some ^-optimal control for (3.1) with the cost (2.6) for which the solution to (3.1) is unique in the weak sense. This assertion would be true in a continuous control case. Since any singular control can be approximated by a continuous control, we deduce that this is true. The above properties mean that our model satisfies the assumptions that were made in [16] to get convergence. We use kh to denote the optimal control for (3.23). In [16], Theorem 5.3 says that kh converges weakly to some admissible control k(-) and the approximate Markov chain under the control kh converges weakly to a continuous process under control &(•). Moreover, Theorem 5.4 says that J*t(&) = V\X) -> JxCK-)). (4.32) At the same time, Theorem 5.6 tells us that V\X) -+ V{X). (4.33) (4.32) and (4.33) show that Jx(k(-)) — V(X). This means that &(•) is an optimal control. These observations say that our Vh(X) and kh(X) are approximations of the real V(X) and the optimal k(X). Chapter 5 Test Results In (2.1) and (2.5), we have parameters a i , a 2 , 6n, 612, &215 2^25 c, fi and a. We pick some parameters to test the numerical methods suggested in the previous chapters and show the results. We do not consider the economic meaning of the parameters in this chapter. As is stated in chapter 4, there are two classical numerical methods: iteration in policy space (IPS) and iteration in state space (ISS). Under the category of state space iteration, we have Jacobi and Gauss-Seidel methods. The numerical tests we did using the above mentioned methods show that: 1. The running time of IPS and ISS are the same, although the iteration number for IPS is much smaller than that for ISS. The explanation is that at each iteration step of IPS, there is a big linear system to be solved which consumes the most computer time. 2. The Gauss-Seidel method is approximately twice as fast as Jacobi method; this is a well known numerical result. Whatever the method we use, one thing is certain: they all lead to the same numerical approximation. So in the following discussion, we are not going to mention the method we chose. Our numerical results show that the control solution divides the whole area into three regions: i) a so called "no-control" region. Actually when the process reaches the boundary of this region (we call it "free boundary"), the control will push the process a little bit, just enough to keep the process inside the region. Otherwise, the control is not applied, ii) a region where the control pushes the process right into the "no-control region" in one step, and then stays in the "no-control region" thereafter; iii) a region where the control pushes the process left into the "no-control region" in one step, and then stays in the "no-control region" thereafter. 32 Chapter 5. Test Results 33 Another point we have to keep in mind when we do numerical tests is the restriction (3.18) on the covariance matrix of the stochastic process. Throughout this chapter, the coefficients are such that (3.18) is satisfied. 5.1 Test 1: uncontrolled problem. We got the exact cost function for the uncontrolled process and compared it with the results from the Markov chain approximation. In the uncontrolled problem, c = 0 and a = 0 . For simplicity, we let &12 = 0 and p, = 0. So we are solving the following problem dxi(t) = (a-i + bnx^t^dt + aidWtt) (5.1) dx2(t) = (a2 + b2ix1(t) + b22X2(t))dt + a2dW(t) with initial values (2.2) and cost function Jx(0) = Ex | ^ ° ° i a ^ e - " * * } . (5.2) The notation in (5.1) and (5.2) is the same as in (2.1) and (2.5). The way to solve (5.1) and (5.2) analytically is illustrated in Appendix C. We use ai = 0.05, bn = - 1 , a2 = 0.1, b2X = - 1 , b22 = - 3 , ax = (0.03,0.01), a2 = (-0.01,0.02) and p = 1 with grid size h = 0.005 and h = 0.01 over a region [-0.2,0.3] x [-0.2,0.3]. That is, An = -0 .2 , A12 = -0 .2 , A21 = 0.3 and A22 = 0.3 (cf. Figure 3.1). This region is shown in Figure 5.1. The results are shown in Figure 5.2, Figure 5.3 and Figure 5.4. When reading Figure 5.3 and Figure 5.4, please refer to Figure 5.1 for the names of the different lines in the region. Some observations can be made from these figures: 1. The Markov chain approximation gives a very close estimate of the real solution. 2. As the mesh size h gets finer, the error in the approximation gets smaller. Figure 5.3 and 5.4 show that the error with h = 0.01 is noticible but with h = 0.005 is negligible. This shows the convergence of the approximation. Chapter 5. Test Results 34 top line ( -0.2, 0.3) (0.3, 0.3) vertical ____ middle line right line left line v horizontal middle line (-0.2, -0.2) (0.3, -0.2) bottom line Figure 5 . 1 : The region [ - 0 . 2 , 0 . 3 ] X [ - 0 . 2 , 0 . 3 ] 3 . Errors at the points on the numerical boundary are bigger than those at the interior points. 5.2 Test2: diagonal drift matrix. In this test, we let 612 = h2\ — 0 . This means that the drift in the x\ direction is independent of x2 and vice versa. Since the control is only exerted in the x\ direction, the decision whether to or not to apply the control is not going to influence the process in the x2 direction. Consequently, the control is not going to influence the cost arising from the x2 direction. Hence we expect the region in which the control is active to be perpendicular to the x\ direction. This is exactly what we get from our numerical test; We used a\ = 0 . 0 5 , &n = — 1 , a2 = 0 . 1 , 6 2 2 = - 3 , a i = ( 0 . 0 3 , 0 . 0 1 ) , t r 2 = ( - 0 . 0 1 , 0 . 0 2 ) , ji = 0 . 1 , a = 0 . 0 0 3 and p = 1 on the region [ - 0 . 2 , 0 . 3 ] x [ - 0 . 2 , 0 . 3 ] with h = 0 . 0 0 5 . The control is shown in Figure 5 . 5 . Chapter 5. Test Results 35 x 2 - ° - Z - 0 - 2 X1 exact solution to the uncontrolled problem ^ - 0 . 2 - 0 . 2 X1 numerical solution to the uncontrolled problem (h = 0.01) numerical solution to the uncontrolled problem (h = 0.005) Figure 5.2: Cost function on the region [-0.2,0.3] X [-0.2,0.3] Chapter 5. Test Results Jx(0) and Jx(0) : X = (xlt - 0 . 2 ) Jx(0) and Jx(0) : X = - 0 . 2 ) -0.2 -0.15 -O.I -O.OS O COS 0.1 0.16 0.2 0.26 0.3 h = 0.01 Jx(0) and Jx(0) : X = ( s i , 0 .05 ) -0,2 -0.1S -0.1 -0.05 O O.OS 0.1 0.15 0.2 0.2S 0.3 h = 0.005 Jx(0) and Jx(0) : X = (a: i ,0.05) -0.2 -0.15 -0.1 -0.05 O 0.05 0.1 0.15 0.2 0.25 0.3 -0.2 -0.15 -0.1 -0.05 O 0.05 0,1 0.15 0.2 0.25 0.3 h = 0.01 Jx(0) and Jx(0) : X = (« i , 0 .3 ) -0.2 -0.15 -0.1 —O.OS O 0.05 0.1 0.15 0,2 0.25 0.3 hr'= 0.005 Jx(0) and Jx(0) : X = (xu0.3) 2 -0.1 S -0.1 -0.05 O 0.0S 0.1 0.1S 0.2 0.25 0.3 h = 0.01 h = 0.005 Figure 5.3: Comparison of different grid sizes on horizontal lines Chapter 5. Test Results Jx(0) and j £ ( 0 ) : X = (-0.2, x 2 ) <M°) and j £ ( 0 ) : X = (-0.2, z 2) -0.2 -0.15 —O.I -0.05 O O.OS 0.1 0.15 0.2 0.25 0.3 h = 0.01 Jx(0) and j £ ( 0 ) : X = (0.05, x2) -0.2 -0.15 -0,1 -O.OB O O.OS 0.1 0.15 0.2 0.25 0.3 h = 0.01 Jx(0) and Jx(0) : X = (0.3, x2) -O.Z -0.15 -0.1 -0.05 0 0.05 0.1 0.15 02 0.25 O. .2 -0.16 -0.1 -O.OS O 0.05 0.1 0.15 0.2 0.25 0,3 h = 0.005 Jx(0) and Jx(0) : X = (0.05, x2) -0,2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0,2 0.25 0.3 h = 0.005 Jx(0) and Jx(0) : X = (0.3,x2) -0.2 -0.15 -0.1 -0.05 0 0.05 0,1 0.15 0,2 0.25 O. h = 0.01 h = 0.005 Figure 5.4: Comparison of different grid sizes on vertical lines Chapter 5. Test Results 38 0.31 1 1 1 1 o.25 iiiiiiiiiii ^ ^ ^ ^ 5 l » i l ^ ^ ^ ^ ^ ^ ^ 0.2 IIIIIIIIIII Jill 0.15 IIIIIIIIIII |§§I 0.1 Ijllllllll fijjjjj 0.05 | o IIIIIIIIIII BB -0.05 -0.1 IIIIIIIIIII IB -0.15 IIIIIIIIIII _Q2^ 1 I I ' ' ^0.2 -0.15 -0.1 -0.05 0 0.05 0.1 X1 grey: control to the right black: control to the left Figure 5.5: Control for uncorrelated problem 5.3 T e s t 3 : different grid size. In this test, we use a\ = 0.05, h\\ = —1, 612 = 0.5, a2 = 0.05 , 62i = —1, 622 = — 4, o\ = (0.03,0.01) and a2 = (-0.01,0.02) with p = 1, fj, = 0.1 and a = 0.002. We calculate the control using different grid sizes h = 0.01, h = 0.005, h = 0.001, and compare the free boundaries thus derived. The results are in Figure 5.6 and Figure 5.7. It shows that as h decreases, the numerical solutions seem to converge. (Nb. The free boudaries for h = 0.005 and h = 0.001 are almost identical.) 5.4 Test 4: heuristic discussion of the relationship between the deterministic dif-ferential equations and the control. In this test, we have two cases. The parameters used for the two cases are in Table 5.1 and we used the grid size h = 0.005. The control for the two cases is shown in Figure 5.8. The difference between the two cases is in the coefficients for the drift part. This difference leads to a different control region: the "no-control strip" in case 1 has positive slope while for case 2 it has negative slope. IK HHHHHHHI 0.15 0.2 0.25 0.3 Chapter 5. Test Results 39 0.25 0.2 0.15 0.1 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 XI h = 0.01 0.3 * o.25 l l l l l l l l i i i l l l l l l l l l l l l l l l l i i l j l o . 2 i ! i i i ! ^ 0.15 1 -^0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 x i h = 0.005 -0.2 -0.15 -0.1 -0.05 O 0.05 0.1 0.15 0.2 0.25 0.3 X1 h = 0.001 grey: control to the right black: control to the left Figure 5.6: Comparison of different grid sizes Chapter 5. Test Results 40 - O A S -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 .2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 h = 0.01 kh= 0.001 h = 0.005 kh = 0.001 grey: control to the right black: control to the left Figure 5.7: Comparison of free boundaries for different grid sizes a\ 6n &12 a-2 &21 &22 o-u 012 Cr 2l 0~22 P t1 a case 1 0.05 - 1 0.5 0.23 -2.5 - 4 0.03 0.01 -0.01 0.02 1 0.1 0.002 case 2 0.05 - 1 0.5 -0.09 2.5 - 3 0.03 0.01 -0.01 0.02 1 0.1 0.002 Table 5.1: Parameters for test 4 To understand the reason for this difference, we analyzed the corresponding deterministic linear differential equations: A-rJi\ 0.05 -x1(t) + 0.5x2(t) 0 .23- 2.5 Xi(t) - 4x2(t) ( dx^t) dt dx2(t) (5.3) and dt ( dxx{t) dt dx2(t) dt (5.4) 0.05 - z i ( i ) + 0.5 x2(t) -0.09 + 2.5 x^t) - 3x2(t). (5.3) and (5.4) are autonomous systems where the independent variable t does not appear explicitly. By setting the right hand sides of (5.3) and (5.4) equal to zero, we get the same equilibrium points for (5.3) and (5.4): (0.06,0.02). The X coefficient matrices for (5.3) and Chapter 5. Test Results 41 ^0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 -0 2 - 0 15 -0 1 - 0 05 0 0 05 0 1 0.15 0.2 0.25 0.3 xl xl case 1 case 2 grey: control to the right black: control to the left Figure 5.8: Control for test 4 (5.4) are - 1 0.5 - 1 0.5 -2.5 - 4 2.5 - 3 The eigenvalues for the first matrix are —3.5 and —1.5 with eigenvectors ( — 1,5) and (—1,1) respectively. And the eigenvalues for the second matrix are —3.4 and —0.5 with eigenvectors ( —1,5) and (1,1) respectively. From the theory of ordinary differential equations, we know that the equihbrium points are stable because of the negative eigenvalues. Comparing Figure 5.8 and Figure 5.9, we see that the control pushes in the direction in which the trajectory goes to the equihbrium point slowly, i.e. the no-control regions are transverse to the direction of "slow" attraction. In case 1, this direction is ( — 1,1) and in case 2, it is (1,1), i.e. it is the direction determined by the eigenvector corresponding to the larger eigenvalue. The trajectories of (5.3) and (5.4) in the phase planes are shown in Figure 5.9. This test tries to give some intuition on how the control region relates to the corresponding deterministic differential equations. That is why we call it "heuristic". Up to now, we just picked some numbers for the parameters in the model to test our numerical approximation method. The parameters we have used are for the purpose of test, so Chapter 5. Test Results x2 case 2 Figure 5.9: Phase plane Chapter 5. Test Results 43 they do not involve any. economics. Next we are going to describe our parameter estimation from "real data" and the results. Chapter 6 Data Fitting In (2.1), we need to estimate parameters ai, bn, b\2, a2, 621, b22, c> o"i and o2 to give the model economic meaning. This is the topic for the first part of this chapter. After the estimation, we run the model with the estimated parameters and discuss the results. 6.1 Description of the data. A l l the data are from a data base called C A N S I M , which is produced by Statistics Canada. It is available at the web site http://www.datalib.ubc.ca/finance/cansim. We use the 3-month Canadian Treasury Bi l l rate as the nominal interest rate x\, we take the Bank of Canada's rate as the bank rate k(t), and we calculate the inflation rate x2 from the Canadian Consumer Price Index. The frequency of the data is monthly and the data cover a period of six years from January 1983 to December 1988 as depicted in Figure 6.1, Figure 6.2 and Figure 6.3. 10 20 30 40 SO SO TO Jan83 Jon94 Janes JenBB JanB7 Janes DecBS Figure 6.1: Canadian treasury bill rates 44 Chapter 6. Data Fitting 45 1 1 5 1 0 2 0 3 0 A O S O 6 0 7 0 J a n 8 3 j Q n ( 3 4 J a n 8 5 J a n B B J a n B 7 J a n S S D o c 8 8 Figure 6.2: Consumer price index 1 0 2 0 3 0 A O S O S O 70 J a n 8 3 J o n 8 4 J a n 8 5 J o n S O J a n 8 7 J a n 8 8 D 0 0 8 8 Figure 6.3: The Bank of Canada rates 6.2 E s t i m a t i o n method . Since the data we have are in discrete form, we need to discretize the continuous time process before using an estimation method. A straightforward way is to change all cPs to A ' s in (2.1). So a discrete-time approximation is A z i ( f ) = (a1 + b11xl(t) + b12x2(t))At + a1AW(t) +cAk(t) (6.1) Ax2{t) = (a2 + b21x1(t) + b22x2(t))At + a2AW(t) where AW(t) = \ZA~t(ei,e2)' and the e,-'s {i = 1,2) are independent random variables with standard normal distribution, because W(t) is an 7£ 2-valued standard Brownian motion. Now the question is how good is this approximation. If X(t) is the linear interpolation of the discrete time process X(ti),(i = 1,2,...) that satisfies (6.1), and Xc(t) is the continous process that satisfies (2.1), then we have (cf. Kloeden and Platen [14], p.343) £ s u p \Xc(t) - X(t)\ < I(VA1 t Chapter 6. Data Fitting 46 10 20 30 AO SO SO 70 Jon03 Jan64 JonOS Jan86 Jon87 Jan88 D 0 0 8 8 Figure 6.4: The inflation rates where K is a constant. This shows that X{t) is a strong approximation of order 1/2. And for each smooth function g, we have (cf. Kloeden and Platen [14], p.473) sup \E{g(Xc(t)) - g(X(t))}\ < KgAt t where Kg is a constant depending on g. This shows that X(t) is a weak approximation of order 1. We use the notation xi(i) for £i(2,), Axx(i) for Axi(ti), x2(i) for x2(ti), Ax2(i) for Ax2(ti) and k(i) for k{ti). If CPI(i) stands for the consumer price index at time we calculate x2(i), which is the inflation rate, by x2(i) = (CPI(i) - CPI(i - l))/CPI(i - 1). The inflation rates thus derived are shown in Figure 6.4. We derive Ax-i(i), Ax2(i) and Ak(i) using Axi(i) = x\(i + 1) — Xi(i), Ax2(i) = x2(i + 1) — x2{i) and Ak{i) — k{i) — k(i — 1) respectively. We use a different way of calculating Ak(i) because it is to be a control variable, hence adapted. If we let y3(i) = Axj(i)(j = 1,2), /(t) = Ak(i), CXJ = o,-At(j = 1,2), 0kj = bkjAt(k,j = 1,2) and jkj = 0kjVA~i(k,j = 1,2), then we can write (6.1) as y i ( 0 = « i + Pnxi(i) + Pi2X2(i) + cf(i) + (711,7i2)(£i, £2)' , x (6.2) 2/2(0 = a 2 + /32ia;i(0 + /322a;2(0 + (72i,722)(£i,£2)/-(6.2) is a linear statistics model. We make different assumptions on Pi2 and c and compare the results. Altogether we have four cases as shown in Table 6.1. We use the method of moments for estimation. We are going to use case 3 for the explanation of our method, but the idea is exactly the same for the other cases. Since c = 1 in case 3, we Chapter 6. Data Fitting 47 012 c case 1 fix 0 1 2 = 0 fix c = 1 case 2 fix 0 1 2 = 0 c is to be estimated case 3 012 is to be estimated fix c = 1 case 4 /?12 is to be estimated c is to be estimated Table 6.1: Four cases (6.3) move f(i) to the other side of the equation and let y[(i) = yi(i) — f(i), then we have 2/i(0 = a i + 0 n * i ( O + (312x2(i) + ( 7 n , 7 i 2 ) (£ i , £ 2 ) ' y2(i) = a2 + 0 2 i * i ( O + 022*2(O + (72i,722)(£i ,£2)'. Since £\ and s2 are mutually independent standard normal variables, and are independent of Xi and x2, we have the following equations if we let £i = y[ — a\ — fl\\X\ — (i\2x2 and 6 = 2/2 - OL2 - 0 2 i * i - fl22x2: E^ = 0 E[(Zi ~ E^)(Xl - EXl)] = £ ( 6 * i ) = 0 E[({i - E^)(x2 - Ex2)} = E(^x2) = 0 E£2 = 0 E[(h ~ E^)(x! - EXl)] = £ ( 6 * i ) = 0 E[({2 ~ Et2)(x2 - Ex2)} = E(Z2x2) = 0 ^ i 2 = 7?i + 7 i 2 2 £ f l 6 = 711721 + 712722 E£,2 = 721 + 7 2 2 -If we use "average" for the expectation, then we have N X ) (y i (0 - « i - 0 n * i ( O - Pi2x2(i)) = o N N [ (» i (0 - "1 - 011*l(O - 012*2(0) *!(*)] = 0 N 1 i=l N jj E M(i) -Cti- )9l lXl(*) " 012*2(0) ^2(0] = 0 i=l N 1 E(^ 2(*) ~ 0 1 2 ~ 021*1(0 - 022*2(0) = 0 4 = 1 (6.4) (6.5) (6.6) (6.7) Chapter 6. Data Fitting 48 1 N J £ [(2/2(0 - "2 - 0 2 1 * l ( O - 022* 2 (O)* i (O] = 0 (6-8) i=i 1 ^ X £ [(2/2(0 - « 2 - feariC*) - 0 2 2 * 2 ( O ) * 2 ( O ] = 0 (6.9) i=l 1 " ]yT E (2 / i (0 - « i - / ? i i* i (0 - & 2 Z 2 ( 0 ) 2 = 7 u + 712 (6.10) i=l 1 N Jf £ [(»l(0 - ^ 1 - 0 1 1 * l ( O - 0 1 2 * 2 ( 0 ) 8 = 1 1 N (2/2(0 - « 2 - 0 2 1 * l ( O - 022*2 (0 ) ] = 711721 + 712722 (6.11) N E(2/2(0 - « 2 - 0 2 1 * 1 ( 0 - 0 2 2 * 2 ( O )2 = 7221 + 7222- (6-12) i=l We solve (6.4)-(6.6) for a u 0 n and / J 1 2 and (6.7)-(6.9) for a 2 , 0 2 i and / J 2 2 . (6.10)-(6.12) give the elements in the covariance matrix of £1 and £2- We assume that 712 = 7 2 i , then we can find the solutions oi~jkj{k,j = 1,2). This assumption of symmetry does not influence our result, because in practice, the information we really want is the covariance matrix of £1 and £ 2 , that is, the information about how £1 and £ 2 correlate with each other, rather than the information about each process itself. Also in the Markov chain approximation in Chapter 3, all we need is the covariances ay's = 1,2). The ay's = 1,2) never show up explicitly in the method. The final step is to change the ctj (j = 1,2), fikj (k,j = 1,2) and jkj (k,j = 1,2) back to aj (j = 1,2), bkj (k,j = 1,2) and cr/y (k,j = 1,2). The results are in the next section. 6.3 Estimation Results. The results are shown in Table 6.2. Note that we also calculate the equilibrium points, which are derived in the same way as in Section 5.4. While the estimation method we have used is fairly basic, the results we obtain do display some meaningful characteristics. First, we have almost the same equiUbrium points in all four cases. Second, the values of the cr/y's ((k,j =1 ,2 ) give us a diagonally dominant covariance matrix so that the restriction on our numerical method (3.18) is satisfied. Third, the sign of frfcj's (k,j = 1,2) shows a certain pattern: &n, 6 2 i and 6 2 2 are negative while &12 positive. This Chapter 6. Data Fitting 49 case 1 case 2 case 3 case 4 ai 0.171 0.116 0.149 0.084 611 -1.805 -1.215 -1.705 -1.058 bi2 0 0 0.307 0.446 c 1 0.354 1 0.341 &2 0.971 0.971 0.971 0.971 hi -4.151 -4.151 -4.151 -4.151 h2 -13.974 -13.974 -13.974 -13.974 O i l 0.0199 0.0156 0.0199 0.0152 012 = 0-2l 0.0005 0.0004 0.0005 . 0.0004 022 0.0830 0.0830 0.0830 0.0830 equilibrium point (0.095,0.041) (0.095,0.041) (0.095,0.041) (0.096,0.041) Table 6.2: Estimation results pattern makes the equilibrium point of the deterministic differential equation a stable one. The reason is the following. According to the theory of ordinary differential equations, the stability of the equihbrium point is determined by the following 2 eigenvalues: (hi + 622) ± \/(bu + 6 2 2) 2 - 4(fen62  ~ buhl) 2 The part under the square root is always less than |6n + 62212 because 6n6 2 2 — 612621 > 0 (cf. point three above). If the part under the square root is nonnegative, then the eigenvalues are real, negative. If the part under the square root is negative, then we have two complex eigenvalues with negative real part. Both cases give us a stable equilibrium point. A stable equihbrium point means that over the long run, the interest rate and inflation rate tend to return to this point if perturbed away from it. This is one of the properties that a stable economy should possess. But when we did estimation using the data on a longer time interval, for example, from July 1982 to December 1988, or even a longer period, from October 1978 to October 1994, the sign pattern we mentioned in point three above does not hold, although the equihbrium points were still close and the covariance matrices were diagonally dominant. The major change was that the sign of 621 turned positive. 621 measures the effect of the interest rate on the inflation rate. If 621 is positive, then high interest rate will push up the inflation rate. If 621 is negative, Chapter 6. Data Fitting 50 then high interest rate will drive down the inflation rate. Practically, the Bank of Canada has operated under the assumption that 621 is negative for the past fifteen years. But it is not clear at this stage what causes the sign difference of 62i in our estimation. It might be that the data we use are not accurate or appropriate. It might be that we need a more advanced estimation method or it might be that we need to modify the original continuous stochastic model. The research is ongoing. Right now, we consider the situation when 621 is negative. 6.4 Control results and interpretation. We use the grid size h = 0.005 for all the numerical calculation in this section. Figure 6.5 shows the control in the four cases with p = 0 .1 , a = 0.002, p = 0 .1. The control regions look very similar to each other, so we use the parameters of case 4 to do the remainder of the tests. Figure 6.6 compares the control for case 4 with different p's: 0 .1 , 0.05 and 0.02 and [i = 0 .1 , a = 0.002. As p gets smaller, more control shows up, but the difference is not significant. Figure 6.7 compares the control for case 4 with different /i 's: 0, 0.1 and 0.5 and a = 0.002, p = 0.05. The graph for fi = 0 is not very clear because we only get one side of the control. To see the whole picture, we run our program on a bigger region [—0.2,0.55] X [—0.45,0.3] . The result is shown in Figure 6.8. When we compare these graphs, we can see that as fi gets smaller, the "no-control strip" not only gets wider but also turns clockwise. A smaller p means that we care less about what x\ is, so we use less control to move X\. We would rather use control to get at x2-Figure 6.9 compares the control for case 4 with fi = 0 .1 , p = 0.05 and a taking values 0.002 and 0.005. As we can see, as a getting bigger, less control shows up because larger a raises the relative cost of the control. As we said before, we use this model to study the dynamics of the nominal interest rate and the inflation rate, with a central bank trying to control the inflation rate by influencing the nominal interest rate. From this point of view, the model with parameters a = 0.002, p, — 0.1 Chapter 6. Data Fitting 51 025 0 2 .6 1 " I 0 0 5 0.2 0.25 0.3 case 1 case 2 s 0 1 s i 005 is 1 0 -0 05 -0 1 -0 15 -0.2^ S 005 0.2 0.25 0.3 case 3 case 4 grey: control to the right black: control to the left Figure 6.5: Comparison of the four cases (p = 0.1) and p = 0.05 give us an appealing result. Please refer to Figure 6.6 (case 4: p = 0.05). Since usually we have nonnegative interest rates and inflation rates, let us look at the graph in the region [0,0.3] X [0,0.3]. What this picture tell us is the following: • When the interest rate is relatively low, the inflation tends to go up. So the bank tries to stem this tendency by raising the short term interest rate. This corresponds to the grey area in the graph. • When the interest rate is relatively high and the inflation is in a reasonable range, the bank wants to stabilize the economy by lowering the interest rate. This corresponds to Chapter 6. Data Fitting 0.25 O.i. -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 interest rates case 4: p = 0.1 0 .3 j : 0 . 2 5 tllltl l l l l l l i i l i.i i i i l l ! interest rates case 4: p = 0.05 0.3 0.25 0.2 -0.15 jllJIf llf?' ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ M _ Q 2 Li I I I I I i \nifi 11 II T ftfflffftiifrfrfrffiiftftiiihiii -0.2 -0.15 -0.1 -O.05 O 0.05 0.1 0.2 0.25 0.3 interest rates case 4: p = 0.02 grey: control to the right black: control to the left Figure 6.6: Comparison of different p's Chapter 6. Data Fitting 53 case 4: p = 0 case 4: p = 0.1 0.3 0.25 0 2 0.15 S .5 0 0 5 -O 15 -0 2 -0.2 -0.15 -0.1 -0.05 O case 4: p, = 0.5 grey: control to the right black: control to the left Figure 6.7: Comparison of different /x's Chapter 6. Data Fitting 54 I 0 1 | 0.0! ~°-0.2 -0.15 -0.1 -O.OS 0 0.05 0.1 0.15 0.2 0.25 0.3 mterest rales case 4: p = 0 on region [-0.2,0.3] X [-0.2,0.3] ^mimmm mmm^ ;; wm^^M 0.2 -0.3 h -0.4 interest rates case 4: \i = 0 on region [-0.2,0.55] x [-0.45,0.3] the smaller box is region [-0.2,0.3] x [ -0 .2 ,0 .3 ] grey: control to the right black: control to the left Figure 6.8: Control on a bigger region Chapter 6. Data Fitting 55 s e § 0 05 - 0 2 - 0 1 5 -0.1 -0.05 0 0.15 0.2 0.25 0.3 case 4: a = 0.002 -0 15 -0.2' -0.2 -0.15 -0.1 -0.05 0 case 4: a = 0.005 grey: control to the right black: control to the left Figure 6.9: Comparison of different a's the black area in the graph. • There is a region where the bank tends to do nothing (corresponding to the white area in the graph), either when the interest rate and the inflation rate are in a reasonable range, (which means that the economy is healthy at that moment, so no adjustment is needed) or when both the interest rate and the inflation rate are high (this is a very complicated economic phenomena and there are different explanations, but a common view is that under such circumstance the interest rates are high enough to drive down the inflation rates.) Another comment we can make is about p. As we can see, the number we take for p is small. What this says is that, although the bank wants to keep both the interest rate and the inflation rate in a certain range, it would be willing to sacrifice the interest rate to control the inflation rate. C h a p t e r 7 C o n c l u s i o n In this thesis, we used the Markov chain approximation method to solve a central bank control problem. In our model, the central bank changes the bank rate to affect the interest rate aiming to control the inflation rate. In terms of mathematical language, our model, mainly from the work of Chiarolla and Haussmann [5], is a two-dimensional stochastic singular control problem with infinite horizon. We constructed a Markov chain to approximate the continuous process and established the corresponding discretized control problem by following Kushner and Martins [16] and Kushner and Dupuis [15]. We derived a dynamic programming equation from the discretized control problem and solved the dynamic programming equation using two iteration methods, namely, the iteration in policy space and the iteration in value space. The convergence of the iteration methods and the convergence of the Markov chain approximation were discussed heuristically. Then we tested our methods on some examples for which we know the results. Those test results confirmed our convergence discussion and showed that the optimal control divides the state space into a "no-control" region and a "control-active" region. The "no-control" region is where the control is not applied when the process stays inside the region. But when the process reaches the "free boundary", the control will push the process a little bit to keep it inside the region. So we actually abuse the word "no-control" here. The "control-active" region is where the control pushes the process into the "no-control" region instantly and the process stays in the "no-control" region thereafter. The "control-active" region has two subregions: one is where the bank rate is raised to increase the interest rate and the other is where the bank rate is reduced to lower the interest rate. The "no-control" region is "hourglass" shaped and squeezed 56 Chapter 7. Conclusion 57 between the two "control-active" subregions. Also, we showed heuristically how the slope of the "no-control hourglass" relates with the corresponding deterministic differential equations. At the end, we used some "real" data from Statistics Canada to estimate the parameters in our model and run our program with the estimated parameters. Then we interpreted the control results thus obtained in an economic context. As we have noted, different time series of the data gave us different patterns of the param-eters and hence different control results. At present, we are not clear about the reason for this. This gives us some hints for future work. For example, we might need to look for other data sources for the interest rates and the inflation rates. Also, we might need to investigate a better estimation method. Another possibility is that the model itself does not reflect the dynamics of the real financial world very well. As is mentioned at the end of Chiarolla and Haussmann [5], maybe we should take the influence of the foreign currency exchange rates into consideration. Generally speaking, our current model is a starting point for studying the central bank control problems, and the numerical solutions gives us a vivid illustration of how realistic or not the model is. Appendix A Derivation of the Transition Probabilities for the Diffusion Case We mentioned in Section 3.2 that we constructed the transition functions for the diffusion case by using a finite difference approximation to the derivatives in the uncontrolled problem's H J B equation (3.4). In this appendix, we are going to show how to derive (3.9)-(3.16) from the difference scheme (3.5), (3.6), (3.7) and (3.8). The notation in this appendix is the same as in Section 3.2. Suppose we already have a "locally consistent" Markov chain with transition probabil-ities ph(X,Y) (X,Y G Gh) and time step Ath(X), that approximates the uncontrolled process, then by the same argument as in the first part of Section 4.1, the value function Vh(X) at point X for this Markov chain satisfies the following dynamic programming equation: Vh(X) = e- pAth^J2ph(X,Y)Vh(Y) + ^(fixl + x22)Ath(X) (A.l) Y 2 Now our objective is to find ph(X, Y) ( X , Y e Gh) and Ath(X). We substitute (3.5), (3.6), (3.7) and (3.8) into (3.4), collect terms, let the term involving V(X) and the other terms stay on different sides of the equation, and note that a\2 = « 2 i and that (3.5) implies B,(x)d-^1 -~ B + ( x ) n * + ^ ) - y f f l _ y ( * ) - n * - « * ) oxi 1 h % h where Bf(X) = max[5,-(X),0] and B~(X) = max[-Bi(X),0]. We can write (3.7) and (3.8) in a similar way. Then we have p - K i - (a+ + ai"2) + a22 + (Bt + Bi)h + (B+ + B2)h + ph2]V(Xl,x2) = / M l a l l - 2a12 _ la12 +  h B l ]Vi xl +  h, x2) + 7^ [|«u - | « i 2 - 2 a i2 + hB^]V(x1 - h,x2) 58 Appendix A. Derivation of the Transition Probabilities for the Diffusion Case 59 + ^[\a22 ~ \at2 ~ \ a \ 2 + hB%]V(x1,X2 + h) + p-[f«22 - | a i 2 - \o-i2 + hB2]V(x1,x2 - h) + / ? 2 a i 2 ^ ( x i + h'->x2 + h) + ^ 2 a t 2 v ( x i - h,x2- h) + Wlai2V(xi + h, x2 - h) + ^\a^2V(xi - h,x2 + h) + \(px\ + x\). (A.2) If we use Qh(X) as defined in (3.17), then the coefficient of V(xi,x2) becomes ^ m + ^ ) = o ^ ) ( 1 + _ ^ L ) . ( A , , If we divided both sides of (A.2) by (A.3) and note that when h is small, ph2/Qh(X) is small, we can say that Also, for a real number r, we have r + + r~ = \r\. Then (A.2) is approximately equivalent to V(Xl,x2) = V(X) = e s7* { QJ^X) [ | « H - \ a l 2 ~ \ a \ 2 +  h B l ]V(X1 + h> Xl) + Q ^ X ) f 2 a n - | a j " 2 _ | a i 2 + fe5f]V(a:i - /i,z 2) + QH^f)[2 a22 - 2 a12 - \ a \ 2 + hB$]V(x1,X2 + h) + ^ X ) ^ 2 2 ~ § a i 2 ~ ^ a i 2 + k B 2 ~ ] V ( a : i , a : 2 - h) + Qf^x)kai2V(xi + h,x2 + h)+ Q7^ X ) | a i " 2 F ( x i -h,x2- h) + Qiqx)kai2V(x~L + M 2 - h) + g^\ai2V(xi - h,x2 + h)} + \(lixl + xl)Q$rje=&. (AA) If we compare the coefficients of V's in (A.4) with those in ( A . l ) , we get (3.9)-(3.16). Note that the difference scheme has been used suggestively for the construction of transition proba-bilities. The important thing of the derived transition probabilities is the "local consistency", which is shown in the next appendix. Appendix B Verification of Consistency Conditions In Section 3.2 we mentioned that our approximate Markov chain satisfies the local consistency conditions (3.20) and (3.21). In this section, we are going to show how (3.20) and (3.21) are derived. A l l the notation in this section is the same as in Section 3.2. First, we show (3.20) is true. According to Figure 3.4, there are eight values for A£^, namely, (h, 0),(-/i, 0),(0, h),(0,-h), (h, h),(—h, -h),(-h, h),(h, -h). So we have EHX,n^hn + + + V h / \ h I ( - h \ h , ph((x1,x2),(xi + h,x2)) + ph((x1,x2),(x1,x2 + h)) + \ ° / ph((x1,x2),(x1 - h,x2)) -h ph((x1,x2),(x1,x2 - h)) -h ph((x1,x2), (xi + h, x2 + h)) + | | ph((x!,x2), (Xl - h,x2- h)) h Ph((x1,x2), (x-i - h,x2 + h)) + i k ^ ph((x1,x2), (Xl + h,x2- h)) If we put (3.9) to (3.16), (3.17) and (3.19) into the above equation, we get (3.20). Next we show that (3.21) is true. Because [ A ^ - ^ , n A ^ ] [ A ^ - ^ i n A ^ ] ' (Atf , i - *&,„At f f l ) a (Atf f l - EXtnAtitl)(Aen,2 - Ex<nA^<2) (A&i - ^ , n A e ^ ) ( A ^ , 2 - 2& i B Atf | 2 ) (Af£ i 2 - 2& i B At f i 2 )2 where A £ ^ ( i = 1,2) denotes the component of the vector A£^. So we are going to prove 60 Appendix B. Verification of Consistency Conditions 61 (3.21) element by element. A ^ j assumes only the values 0, h and — h, and we have ph(A^1 = h) = ph((x1,x2),(x1 + h,x2)) + p\( + Ph(( / ( A ^ j l = - / 0 = p \ + Ph( + Ph{ / ( A ^ 1 = 0) = ph{ + Ph{ x\,x2),(x1 + h,x2 + h)) X!,x2), (xi + h,x2 - h)) xi,x2),(xi - h,x2)) xi,x2),(xi - h,x2 - h)) x\,x2),(xi -h,x2 + h)) xi,x2),(x-i,x2 + h)) x\,x2),(xi,x2 - h)), so E(A^ - £ |Xi ) 2 = E(AtiAf - (EA&Ay = 02ph(A^tl = 0) + h2ph(Atitl = h) + (-h) V ( A # i X = -h) -{Bx(X)At\X)f = anAth(X) + {B^X^At^X^ - Bl(X)(Ath(X))2 . Similarly, we get And o(At»(X)) E(A&2 - EXtnA&2)2 = a22Ath(X) + \B2(X)\At\X)h - B2(X)(Ath(X))2 . y v ' E(A^A - EXtnAtitl)(Aen,2 - EhXtnAen,2) E(AU,ihA^2) - (EA&tl)(EA&2) (h • 0)ph((Xl, x2), (x a + h, x2)) + (-h • 0)ph((xi,x2), (xi - h, x2)) + (0 • h)ph((xX, X2), ( S i , X2 + h)) + (0 • ( - h ) ) ^ ^ Z 2 ) , 532 - h)) +(h • h)ph((x1,x2), (xi + h,x2 + h)) + ((-h) • (-h))ph((Xl,x2), (xx - h , x 2 - h)) +((-h) • h)ph((x1,x2),(xi - h,x2 + h)) + (h • (-h))ph((x1,x2),(xt + h,x2- h)) -(B!(X)At\X))(B2(X)Ath(X)) Appendix B. Verification of Consistency Conditions 62 = a12Ath(X) + (-B1(X)B2(X)(Ath(X))2) = a21Ath(X) + o(Ath(X)). Here we use the fact that a\2 = a2\. These equations have shown that (3.21) is true. Appendix C Cost Function of the Uncontrolled Problem: the Exact Solution The first test we did in Chapter 5 was to compare the exact solution and the numerical solution to the problem defined by (5.1) and (5.2). In this section, we are going to show how to solve (5.1) and (5.2) analytically. The notation here is the same as in Chapter 2. If we switch expectation and integral in (5.2), by Fubini's theorem, it becomes / • O O 1 Jx(0) = / -ExixKt^e-^dt. ( C l ) Now the problem of looking for the cost function is turned into looking for Ex[x\{t)). By Ito 's lemma, d(x\(t)) = 2x1(t)(a1 + bux^tydt + axldt + 2x1(t)a1dW(t) (C.2) where a n is as defined in (2.9). If we write (C.2) in integral form, then we have x\(t) = z 2 (0) + [ (2a1x1(t) + 2buxl(t) + a n ) d f + 2 f Xl(t)dW(t). (C.3) Jo Jo If we take expectation on both sides, then we have Exx\{t) = Exxl(0) + Ex f\2a1x1(t) + 2buxl{t) + an)dt + 2EX f ari(t)dW(<). (C.4) Jo Jo Since Ex\ < oo because xi(t) with xi(0) = x\ is a solution to (5.1) (cf. [12], p.289), we have: EX f' Xl(t)dW{t) = 0 Jo because the stochastic integral is a martingale (cf. [12],pp.137-139). And if we switch the expectation and integral in (C.4), then we can write (C.4) as Exx\(t) = Exxl(0) + f Ex(2alXl(t) + 2bnxl(t) + au)dt, Jo 63 Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution 64 or dEx(xl(t)) = (2a1Ex(x1(t)) + 2b11Ex(x2(t)) + au)dt. (C.5) Similarly, dEx{x\(t)) = (2a2Ex(x2(t)) + 2b21Ex(x1(t)x2(t)) + 2b22Ex(x22(t))dt dEx(xi(t)x2(t)) = (a2Ex(xi(t)) + a1Ex(x2(t)) + ( 6 U + b^Exix^x^t)) +b21Ex(x2(t)) + \{a21 + a12))dt (C.6) (C.7) dEx(xi(t)) = (fll + b11Ex(x1(t)))dt (C.8) dEx(x2(t)) = (a2 + b2iEx(Xl(t)) + b22Ex(x2(t)))dt (C.9) (C.5) through (C.9) is a hnear system of differential equations for Ex(xi(t)), Ex(x2(t)), Ex(x\{t)), Ex(x2(t)) and Ex{x\(t)x2(t)) with initial condition Ex(xi(0)) = X l Ex(x2(Q)) = x2 Ex(xlm = xl Ex(xKO)) = x\ £ x ( x i ( 0 ) z 2 ( 0 ) ) = XxX2 ( C I O ) Starting from (C.8), we can solve (C.8), (C.9), (C.5), (C.7), (C.6) one by one for Ex(xx{t)), Ex(x2(t)), Ex(xl(t)), Ex(Xl(t)x2(t)), and Ex(x2(t)). The result is Ex(x1(t)) = (x1 + ^ ) e b ^ - ^ -On " i i C12 JP / / , u 021C11 b t (a2-b2\C\2 b2\C\\ a2-b2\C\2 Ex{x2{t)) = 7 - e + ( 7 7 7 - + x2) + On - °22 O22 O n - 02  —O22 C21 C22 C23 Ex(x2(t)) = e*»* + {x\ +  a u - J a i ° 1 2 + ^ i ) e2b»t + " n "J j*" 1* 7" 2 o n on - 2 o n ' n C31 C32 C*33 Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution Ex(Xl(t)x2(t)) = J^- eb»< + eb^ + - 5 " e 2 ^ C 4 1 C42 C 4 3 „ O i l + t>22 022 O i l On - 0 2 2 -(On + 0 2 2 ) where = 1(^ 21 + O i 2 ) - a 2 C i 2 + O1C23 + O21C33 D12 = c ^ C i i + 01C21 + O21C31 D13 = «lC"22 D14 = 021C32 and Ex(4(t)) = - r - ^ r - e 6 " ' + ^ + e 2 ^ O i l — ^°22^ C^iS - 1 1 7 2 2-• C51 C52 C 5 3 1 (ElL 1 °22 + '^?3 _ ^ 2 4 ^25 2x 2622< + + OJ, J. + J. OJ. OJ. u u ^ x 2 ) e ' 26 2 2 2622 - O n 022 2611 — 2622 6 n + D ^ e(6n+i»22)< + D 2 1 pu ~ 6 2 2 y ^-2622, C 5 5 C 5 6 -*22 where -D21 = «22 + 2a 2 C 2 3 + 2621C45 Z>22 = 2 a 2 C 2 i + 26 2 iC , 4 i D23 = 2 a 2 C 2 2 + 2 6 i C 4 2 D24 = 2621C43 D25 = , 2621C44. Note that we solve the differential equations when b\\ 7^  622 and 6 n ^ 26 2 2 because this is case for our test problem. Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution 66 Wi th the solution of Ex{x2{t)), we then find the cost Jx(0) as in ( C l ) . The improper integral is well defined because &n and 622 are negative and p is negative so that Ex(X2(t))e~pt is exponentially decaying. We have j /QN _ lr g 5 1 C$2 C53 C54 C55 Cjtf X K ' 2K-b11 + p -b22 + P -2hi + P -26 2 2 + P -bu - 62 2 + p p ' Bibliography [1] Akian, M . , Menaldi, J .L. , and Sulem, A . (1996). On an Investment-Consumption Model with Transaction Costs, working paper. [2] Bancora-Imbert, M . C . , Chow, P.L. , and Menaldi, J .L . (1988). On the Numerical Approxi-mation of an Optimal Correction Problem. S I A M J . Sci. Stat. Comput. 9, No.6, 970-991. [3] Binhammer, H . H . (1982). Money, Banking and the Canadian Financial System, second edition. Methuen, Toronto. [4] Chan, K . C . , Karolyi, G .A. , Longstaff, F . A . , and Sanders, A . B . (1992). An Empirical Com-parison of Alternative Models of the Short-Term Interest Rate. J . Finance XLVII, No.3, 1209-1227. [5] Chiarolla, M . B . , and Haussmann, U . G . (1996). Optimal Control of Inflation: a Central Bank Problem. Working paper. [6] Fitzpatrick, B . G . , and Fleming, W . H . (1991). Numerical Methods for an Oprimal Investment-Consumption Model. Math. Oper. Res. 16, No.4, 823-841. [7] Fleming, W . H . , Zariphopoulou, T. (1991). An Optimal Investment/Consumption Model with Borrowing. Math. Oper. Res. 16, No.4, 802-822. [8] Hansen, L .P . (1982). Large Sample Properties of Generalized Method of Moments Estima-tors. Econometrica 50, No.4, 1029-1054. [9] Haussmann, U . G . , and Suo, W . (1994). Existence of Singular Optimal Control Laws for Stochastic Differential Equations. Stochastics and Stochastics Rep. 48, 249-272. [10] Haussmann, U . G . and Suo, W . (1995). Singular Optimal Stochastic Controls I: Existence. S I A M J . Control and Optimization 33, No.3, 916-936. [11] Karatzas, I., Lehoczky, J.P., Sethi, S.P., and Shreve, S.E. (1986). Explicit Solution of a General Consumption/Investment Problem. Math. Oper. Res. 11, No. 2, 261-298. [12] Karatzas, I., and Shreve, S.E. (1991). Brownian Motion and Stochastic Calculus, second edition. Springer, New York. [13] K i m , J .T. (1993). Estimation of CEV-Type Interest Rate Process. Working paper. [14] Kloeden, P.E. , and Platen, E . (1992). Numerical Solution of Stochastic Differential Equa-tions. Springer, New York. 67 Bibliography 68 [15] Kushner, H . J . , and Dupuis, P .G. (1992). Numercial Methods for Stochastic Control Prob-lems in Continuous Time. Springer, New York. [16] Kushner, H.J . , and Martins, L . F . (1991). Numerical Methods for Stochastic Singular Con-trol Problems. S IAM J. Control and Optimization 29, No.6, 1443-1475. [17] Marsh, T . A . , and Rosenfeld, E.R. (1983). Stochastic Processes for Interest Rates and Equi-librium Bond Prices. J . Finance XXXVI I I , No.2, 635-645. [18] Noble, B . , and Daniel, J .W. (1977). Applied Linear Algebra, second edition. Prentice-Hall, N . J . [19] Parkin, M . , and Bade, R. (1986). Modern Macroeconomics, second edition. Printice-Hall Canada. [20] Shiller, R . J . (1978). Can the Fed Control Real Interest Rates ? from Rational Expectations and Economic Policy edited by Fisher, S., U . of Chicago Press. [21] Tourin, A . , and Zariphopoulou, T. (1994). Numerical Schemes for Investment Models with Singular Transactions. Computational Economics 7, 287-307. [22] Zariphopoulou, T. (1994). Consumption-Investment Models with Constraints. S I A M J . Control and Optimization 32, N o . l , 59-85. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079781/manifest

Comment

Related Items