NUMERICAL SOLUTION OF A CONTROL PROBLEM I N 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 G R A D U A T E STUDIES D E P A R T M E N T OF MATHEMATICS INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A June 1996 © Bing Han, 1996 In presenting degree at this the thesis in partial fulfilment University of British Columbia, freely available for copying of department publication this or of reference thesis by this for his and study. scholarly or thesis for her Department of Mci-the m a f ' The University of British Vancouver, Canada Date DE-6 (2/88) Columbia I further purposes gain the requirements I agree that agree that may representatives. financial permission. of It shall not the an advanced Library shall make permission be granted is understood be for by the that allowed without for it extensive head of my copying or my written 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. ii 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 H J B 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 differential 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 G 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 R 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 h h 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 H J B equation. (Although sometimes the H J B 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 approximation 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 dx (t) x = (ai + b x (t) dx (t) = (a + b ix (t) + b x (t))dt 2 n 2 + b 2X (t))dt + a dW(t) + cdk(t) 1 2 1 1 22 2 1 (2.1) + a dW(t) 2 2 with initial values xi(0 ) = x x (0~) = x - 2 x (2.2) 2 where x\(i) and x (t) represent the nominal interest rate and inflation rate at time t respectively. 2 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) w (t))' is an 7£ -vafued standard 2 2 Brownian motion with w\(t) and w (t) independent of each other, representing uncertainty in 2 the economy. a\ = ( C T U O\ ) and a = (a \ a ) are volatility coefficients. 'a\, 6 , a , b \, b 2 2 2 n 22 are constants. Let / X(t) \ xi(t) , K(t) = x *2(«) 2 1 i \o 2 0 ] a-i + buXi(t) + b x (t) ^ 12 B(X(t)) K 2 a + b xi(t) + b x (t) 2 21 4 22 ( ck(t) 2 / 2 2 22 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: [* B(X(s))ds Jo X(t) = X+ + o f dW(s) Jo + K(t) t for t > 0. The cost function corresponding to the control k is defined as U OO -j POO 1 -(p.x (t) + x (t))e- dt 2 2 1 + aj pt 2 e-^dv (t)j where Ex denotes expectation given the initial condition X = (xi,x ) and Vk{t) = var(& : 2 [0~,2]). (2.5) k 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£ ) and satisfies the 2 following variational inequality min^CV - pV+^(p,xl + xl),mm{a-V ,a + V }^ Xl = 0 Xl (2.7) where CV = (ai + buxx + b x )V 12 2 Xl + (a + b ix 2 2 1 + b x )V 22 2 X2 1 + - dV 2 2 i ai i,j-l Q X . Q 1 (-) 2 X 3 8 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 Markov Chain Approximation 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 n u m e r i c a l b o u n d a r y a n d the n u m e r i c a l g r i d . 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.) W i t h 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) = dx (t) = (a + b x (t) + b x (t)dt 2 (a + b x (t) 1 2 u 21 1 1 + b x (t))dt 12 22 2 2 + o- dW(t) + cdk(t) + dY (t)-dU (t) 1 + a dW(t) + 2 1 1 (3.1) dY (t)-dU (t) 2 2 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) = A i(i 2 = 1,2). This is illustrated i n Figure 3.1. The other symbols i n (3.1) are the 7 Chapter 3. Markov Chain Approximation x 8 2 u ( 11, 22 ) A A ( 21, 22) A 2 A Ul G (A U i A 1 2 ) (A Y 2 1 > A 1 2 ) 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 < (A \ - Au)/h, 0 < I < (A -Ai )/h} without loss of generality, that (A21 — An) and (A — A\ ) are integral multiples of h. Define 2 22 the "extended /i-grid" G% = {(A n (A 22 — A\ )/h 2 + mh,A 12 22 where we assume, 2 2 + Ih) : - 1 < m < (A 21 - A )/h lx + 1, - 1 < / < + 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 x 2 (Au,A ) (A i,A ) 22 2 A A A A A o o o o o o o o A o o 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 o o o o o o o o o o o o (A , (A ,A ) n <-dGh dGi A A A o G.o 22 21 12 Figure 3.2: h-grid Ai ) 2 and G £ W i t h 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\. 10 Chapter 3. Markov Chain Approximation 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: E A& = B{X)At {X) E [A^ x<n + o(At (X)), h Xtn n (3.2) h E Att}' - £* A£*][Atf - Xin iB = aAt (X) + o(At (X)) h h where a is the same as defined in (2.9) and E\ {g,i<n,£ n denotes the conditional expectation given = 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 ) [(V(X) - V(X - h)] /h if Bi(X) < 0 = B(X)(as in (2.4)) and e;(i = 1,2) is the i (3.5) ei where (Bi(X),B (X)) 2 in t h standard basis element n. 2 The difference scheme approximating the second derivatives is the following. If i = j(i,j = 1,2), d V(X) V(X + h) + V(X - h) - 2V(X) dxidxi h 2 ei ei 2 This is a standard central difference scheme for second derivatives. (3.6) Chapter 3. Markov Chain Approximation If i^ j > and d V(X) dx 'Ox' 2 11 = 1,2), 0 2V(X) + V(X + ejh + ejh) + V(X - e h { e h) 3 (3.7) V{X + ejh) + V(X - e^fe) + F ( X + ejh) + F ( X - ejh) 2h To make the above approximation clearer, the right hand part can be rewritten as 2 1 2 f F ( X + ejh + e h) - V(X + ejh) _ V(X + gh) - V(X) I h h F ( X ) - V(X - e,-fe) _ V(X - h) - V(X - ejh - ejh)' •}/h h h 3 ej 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 j < 0 (i, j = 1,2), similarly, we have t d V(X) dxidxj 2V(X) + V(X + e h - ejh) + V(X - ejh + ejh) 2h so o\ V(X + ejh) + V(X - ejh) + V(X + ejh) + V(X - ejh) ' 2h The grid points used is shown in Figure 3.3(b). Note that the diagonal which has points on it 2 { 2 1 + ; 2 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 (cf. Appendix A ) and recalling that a 12 Notation p ((xi, h 12 = a , we get the following transition probabilities. 21 x ), (j/ y )) stands for the transition probability from point (21,2:2) to point 2 (2/1,2/2) and X = 1? 2 (x ,x ). 1 2 time t+ A t (X) h— h —H Figure 3.4: Diffusion step p ((x ,x ),(x 1 + h,x )) = p ((xi,x ),(x - h,x )) = h 1 2 h 2 1 2 p ((x ,x ),(x ,x + h)) p ((x ,x ),(x ,x - h)) h 1 2 1 2 h 1 2 1 + 2 2 f-\^ hB+(X) IQ\X) B;(X) a + f-^ + a 022 h |«2l| + IQ\X) hBt(X)\/Q\X) hB (X) IQ\X) 2 (3.9) (3.10) (3.11) (3.12) (x + h,x + h)) = \a+ /Q (X) (3.13) p ((x x ),( - h , x - h)) = -a+ /Q (X) (3.14) p ((x x ),( -h,x ^ /Q (X) (3.15) ({ ,x ), h P Xl 2 t h u 2 Xl h u 2 Xl 2 2 2 + h)) = h 2 h 2 h 2 Chapter 3. Markov Chain Approximation p ((x ,x ),(x 1 2 h)) = + h,x - h 1 13 2 \a^ IQ (X) h 2 (3.16) where the normahzation factor (to make the p 's sum to 1) is h Q (X) h =a n +a 2 2 - ± ( | a | + | a | ) + fc(l*iPOI + I ^ W I ) . 21 12 (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 an > 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 h 2 At (X) h Q (X) h (3.19) 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 E „A£ h x E (3.20) X,n aAt (X) h 3.3 = B{X)At\X) +o(At (X)) h (3.21) The transition functions: case 2. This is the case when £^ = X 6 Gh and control is exerted. Since the control is "impulsive", we let At (X) h = 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 p ((x ,x ),(x h 1 2 1 h 2 = -h) = l + h,x )\Ak^ = h) = l. 2 p ((x ,x ),(x 1 - h,x )\Akl) 1 2 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£ 15 = h(i = 1, 2) and Af/£- = /i(i = 1 , 2 ) when necessary. { Figure 3.6: Reflection step Let AY// = (AY// , AY£ ) A and AU^ = ( A [ / £ 2 1 ; A f / ^ ) and combine section 3.3 and section 2 3.4, then we can define Y = £ h n AY/ V = 1 h n i=0 3.5 J2 AU? k = h n i=0 ]T Ak/ i=0 The control problem for the chain. We say that a control policy k is admissible if it preserves the Markov property h P{tn+i = Y\g,Ak$,i< = X} = p (X,Y\Ak ). h h n Let At% = At (£%) if n is a pure diffusion step, and set At ^ = 0 otherwise. Define i£ = h Y^i^o Atf. 1 1 Note that {Afc^;} is a discrete time process here and is known at time i£. W i t h initial condition $ = X, the cost function for the approximating chain under control policy k h is { oo i £ oo 5 ~"'"(M(£n,i) + ( t f , ) ) A i * ( t f ) + « £ ^ e ra=0 where (i = 1,2) represents the 2 2 2 n=0 "j p i "|A^| (3.22) J component of ^ and the objective is to minimize j \ Chapter 3. Markov Chain over all the admissible k s h, Approximation 16 and to find the value function V (X) h = mfJ (k ). h x (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 V (X) H in (3.23). A t 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 + x )At (X) 2 (4.1) h 2 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- ^E(V (Cn+i)\U PATH = e- = X) h P M H ^E V X H {^) = e"'*** W £ V {Y)p\X, h Y) (4.2) Y where the sum is over all the states that the chain can reach in one step from X. B y (4.1) and (4.2), the total cost for not applying any control is c - , A t * ( x ) J2P (X, H Y)V\Y) Y where X — + \(pix\ + xl)At\X) L (xi,x ). 2 17 (4.3) 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 + V ( h Xl The optimal cost function V (x\ h + h,x ) (4.4) 2 + h, x ) is not discounted to the present as in the diffusion 2 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 + V (x -h,x ) (4.5) h 1 2 The control direction we choose should be the one that has least cost, so the cost for exerting control is min [V (xi h + h, x ) + ah, V (x h 2 1 - h, x ) + ah} (4.6) 2 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 p (X, Y)V (Y) h | e x )At {X), 2 h 2 [v (xi + h, x ) + ah, V {x - h, x ) + ah } h min = min + \(iix\ + h h 2 - ^ P 0 £ V (x h y p \X, x Y)V {Y) h + h, x ) + ah, V (x h 1 2 1 2 + \(px\ + (4.7) x )At (X), 2 h - h, x ) + ah} . 2 This is called the dynamic programming equation and it is true for all X £ GhWhen 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 V {X) h = E V (£ ) h x 1 = Y p {X,Y)V {Y) h J h (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) Chapter 4. Dynamic Programming and Computation 19 and (4.8) becomes V (X) = V (the nearest point to X on G ) h (4.9) h h This says that when X is on the numerical boundary, we take the value function V (X) the h same as the value function of the nearest grid point. 4.2 Comparison of the dynamic programming equation and H J B equation. By moving V (X) = V (x\, x ) to the other side in equation (4.7), it becomes h min h 2 \e-^ W £ y p (X, ih Y)V (Y) h min V (x h + \{u.x\ + x )At {X) h 2 + h, x ) - V (x , x ) + ah, V (x h 1 h 2 - h 2 1 2 V (X), h - h, x ) - V (x ,x ) + ah} } = 0 h 1 2 1 2 (4.10) The first part in the braces above can be considered as a "diffusion part". If it is divided by At (X), h 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 H J B 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 'V ( h min Xl + h,x ) 2 - V (x!,x ) h 2 + ah V ( h Xl h - h,x ) 2 - V (x x ) h u 2 + ah\ nmmia h + J V^.a-V*) which is very similar to the second part in the H J B equation (2.7) if the superscript h is dropped. Moreover, (4.9) implies that ^ ~\dG n 4.3 h = 0 where n is the unit outward normal. 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). X Q 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\X ) = V (X ). (4.11) We want to substitute (4.11) into (4.7) where V (Xo) shows up in (4.7). According to the h 0 2 h definition of the transition probabilities in Chapter 3, X o has connection with X\, X and X 3 . 2 So there are three equations in (4.7) that have V (X ) in them, that is, V (X\), h h 0 V (X ) h 2 V (X ). h 3 Check V (Xi) h first. The chain can go from X\ to X o and X . B y (4.7), we can write 2 V\X{) = min {e->**l*i) E r ^ , x ) / ( ^ i ^ ) ^ ( ^ ) 0 +6-^ f e (^)(/(x ,Xo)y V (x (x ) + /(x ,x )y 0 1 - h, x ) + ah} hl h2 By (4.11), the underhned part becomes h u 0 (x )) l h2 (p (X X ) f e 2 h hl 2 4 )A*' (Xi), + h, x ) + ah, V (x h i2 f e 1 + £(M4i + with X i = (£1,1, xi ). 2 + p (X ,X ))V (X ). h h 1 2 2 2 and Chapter 4. Dynamic Programming and Computation 21 So we can define p^(X ,X ) 1 Thus V (X ) = p (X ,X )+p (X ,X ). h 2 .(4.13) h 1 2 1 Q is ehminated from the right hand side of V (X ) h in (4.7). h 0 1 Similarly, we can define p~t(X ,X )=p (X ,X ) + p (X ,X ) h 3 2 3 (4.14) h 2 3 0 so that V (Xo) is ehminated from (4.7). h As to V (X ), h 2 the chain might move from X to X o , but it does not stay at X . 2 By (4.7), 2 we can write V\X ) = min 2 { -^t"(X e h^ 2 ) p )V (Y) h {x Y +e-" ( ^p (X ,X )V (X ) Ath x h h 2 0 0 t 1 (- ) 4 15 + l(»x + xl )At (X ), 2 h 2A V {x h 2 + 2 h,x ) + ah,V (X ) + ah} with X = ( x , i , Z 2 , 2 ) . So V (X ) appears in both diffusion and control terms. B y (4.11), the first underlined part in (4.15) becomes 2tl h 2>2 0 h 2 2 0 p\X ,X )V (X ). h 2 0 2 So in a diffusion step, we define ?(X ,X )=p (X ,X ). (4.16) h 2 2 2 0 Similarly, we substitute the second underlined part in (4.15) by (4.11), and define in a control step p*(X ,X \Ak* 2 By (4.16) and (4.17), V (X ) h 0 2 = -h) = l. (4.17) is eliminated from the right hand side of V (X ) h 2 We have to be more careful when we define p (X ,X ) h 2 2 in (4.7). in a diffusion step when X 2 is a corner point as shown in Figure 4.1, (b). As we can see, XQ communicates with X' , XQ and Q Chapter 4. Dynamic Programming and Computation 22 Xg'. By (4.7), V {X ) h 2 = min {e-' A t t ^)EF«^,l«i }/(^,F)^(F) r +e-* V*){jfi{X ,X' )V\X' ) + th 2 Q Q p {X ,X>>)V {X>>) h h 2 +p (X ,Xg')V (X[l')) h (4.18) h 2 + 1(^1,1 + F (a; h a;^, )Ai (X ), /l 2 2 + / i , z , ) + ah, V {X' ) + ah} . h 2 ) 1 2 2 0 By (4.9), we have V\X' ) = Q V (XH) = V (X'f) h = h V (X ). h 2 So the first underUned part in (4.18) can be written as (p (X , X' ) + p\X ,X' + h 2 Q 2 Q) p\X ,XZ'))V (X ). h 2 2 Hence in a diffusion step, we define ?(X ,X ) 2 = p (X ,X^) + p (X ,Xg) h 2 h 2 2 + p (X ,X[;'). h 2 (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 p s. h, (Nb. p = p when boundary points h h 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 x 2 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). W i t h this ordering and the elimination of dG%, we can express (4.7) in vector form. Define a matrix R h = {r\X,Y),X,Y € G) h by r (X,Y) h ( M + l)(iV + 1) states in dG , the dimension of R h h = e~ ( )ph(X,Y). pAth x Because we have is (M + 1)(N + 1) by (M + 1)(N + 1). The Chapter 4. Dynamic Programming and ( A N A 2 2 Computation 23 ) (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 form of R is shown in Figure 4.3. Each block represents one level of the grid. A s we can see. h R h is a nine-diagonal matrix. Then we can write (4.7) as V h = min {R V h h + C \ mm{M V r where the minimum is taken elementwise. V h vector. A n d C h as well. M r h l h = {V (X),X = {\{nx\ + xl)At (X),(x ,x ) h + ah x I , MfV 2 h + ah x 1}} (4.20) £ Gh} is a 1 by ( M + 1)(N + 1) = X e Gh} is a 1 by ( M + l)(N + 1) vector and Mi are defined as the following: U Mi Mr u (M+i)blocks (4.21) L J (M+i)blocks Chapter 4. Dynamic Programming and Computation 24 where 0 U 1 0 L = '•• (4.22) 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 (N+l)x(JV+l) Chapter 4. Dynamic Programming and 4.4 Computation 25 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 iteration. Algorithm for Jacobi Iteration: 1) Choose an initial value function V (X). H 0 2) The (n + l)st iteration is V (X) = min h n +1 { -pAt*(X) e Zy ?(X, + \{px\ Y)V (Y) H N min [ v ^ z i + h, x ) + ah, Vn(x 2 x + z|)A^(X), - h, x ) + ah] } . 2 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,x ) will be on dGf. By (4.9), we use V^(the nearest point to (x\ — h,x ) 2 2 instead. Similar argument can be applied when x\ — 3) If maxx | V /l ra +1 ( X ) - V£(X)\ on Gh) A \. 2 < c, stop and take V£ (X) as the result. Otherwise repeat +1 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 V (X) H N +1 = min { e - ^ P D (ZY<X?(X + \(lixl + xl)At (X),mm h ,Y)V (Y) H + N +1 [v£(xi + h,x ) 2 + ZY>X?(X,Y)VZ(Y)) aft, V ^ z a - h,x ) 2 + 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 3) If maxx \V£ (X) - V£(X)\ +1 < Computation e, 26 stop and take V£ (X) as the result. Otherwise repeat +1 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 V \X) J2?(X, Y)V (Y) = -» W Ath 0 h e 0 + \(px\ Y for + x )At (X) 2 h 2 1 V (X). h 0 2) A t the (n -f l)st iteration, compute the control first. Let k^ +1 = 0, if e-^ WYy(X,Y)V£(Y) + \{iix\ + h x )At (X) 2 h Y < min ( V ^ ( z i + h, x ) + ah, V ^ x i - h, x ) + ah) . 2 2 If the above inequality is not true, then let ki\ +1 V£(xi and, let kl\ = h, when + h, x ) + ah> V„(xi - h, x ) + ah 2 2 = —h otherwise. +1 W i t h the control just defined, solve the linear system e V£+i(X) -^) = I V* ( +1 Xl E y JJ y * (X) + h, x ) + ah 2 Vn (x-i - h,x ) +1 2 + ah )V +i(y)+ + ^ ( j r ) if = 0 if k h n+1 if kl\ +1 = h (4.25) = -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 3) If maxx \V^ (X) - V£(X)\ +1 Computation 27 < e, stop and take V£ (X) +1 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£ (X) in step 2. +1 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 V (X) for each h. The second one is the convergence of h V (X) to the exact value function h 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* = min {R V£ + C , m i n { M V * + ah x 1, M v£ = min {R V£ + C , M V% + ah x 2, MtV* + ah x l } h +1 + ah x 1}} h r h t ^ h r where the minimum is taken componentwise. To explain the concept of iteration matrix, we want to write (4.26) as K + i = R (k )V£ + C {k ) h where R (k ) h n The i t h and C (k ) h n row of R (k ) h and (4.22)) and the n (n = 1,2 ...) h n n (4.27) are composed in the following way. Suppose we consider ( V ^ ) ; . +1 is chosen from the i row of R th element of C (k ) h n (cf. Figure 4.3), M or M h ; is chosen from the element of C r h (cf. (4.21) or the number ah, depending on where the following minimum is achieved: min{ i t h row of R • V +i i t h row of M -V + ah, i t h h h h r t h row of C \ row of Mi • V + ah] h where '•' is a dot product of a row vector and a column vector. So R {k ) h n 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 R (k ) the h n iteration matrix. For each n, R (k ) might be different. For each row of R (k ), h h n n seen, there are three choices, and there are ( M + l)(iV + 1) rows of R {k ). h n there are 3(^+1)^+1) choices for R (k ) (n = 1,2...). This R (k ) h h n n as we have So theoretically 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. B y an argument similar to the one used for the Jacobi method, we have V} +1 where L (k ) h n + +1 U (k )v£ + h C (k ) h n n is the strict lower triangular matrix (not including the main diagonal) of n and U (k ) h = L\k )V^ R (k ), h n is the upper triangular matrix including the main diagonal. Similarly, for the n method of iteration in policy space, we can write (4.25) in vector form vf+i = R (k )v£ + h n+1 c\k ). +1 n+1 The usual convergence argument (such as the one in Kushner and Dupuis [15], pp.154-157) needs the contraction property of R (k ), h n (R (k )y h n We only know that R (k h = 0) = R ^ 0 as i -»• 00. as defined in (4.20) has this contraction property. The h n that is reason is stated below. According to the definition of R , the sum of the absolute values of the elements in a row is h £ \r (X,Y)\ h YeG h = £ e-^W?(X,Y) = -<^W e YeG £ ?(X,Y) YeG h h where Y^YeG P (X, Y) is the probability that X moves to all the state in the state space which h h gives 1. So we have • £ \r (X,Y)\ h YeG h = e- W pAth <1. Chapter 4. Dynamic Programming and This is true for every row, so the H ^ H o o Computation norm of R 29 is less than 1, that is = max £ \r (X,Y)\ h < 1. Yeo X h Let |A| denote the largest magnitude of the eigenvalue of R . Then |A| < h Noble and Daniel [18], p.268). Hence (R y -»• 0 as i -> oo (cf. Noble and Daniel [18], p.378). h For arbitrary R (k ), h R h we can see from the construction of R {k ) n < 1, (cf. Hif^Hoo h n that some of the rows of are replaced either by the corresponding rows of M or the corresponding rows of M ; . From r (4.21) and (4.22), we see that each row in M or Mi has all entries zero except one which is 1, r but that makes ||i2 (A; )|| /l n R (k ) = 1. So the above argument for R h 00 does not apply any more. Also is non-symmetric. Besides that, as we mentioned before, the matrix form of R (k ) h has h n n a huge amount of choice in the minimization. A l l these factors make the study of the properties of R (k ), h n 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 V (X) to the exact value function V(X) h is based on Kushner and Martins [16], Section 5, and Kushner and Dupuis [15], Section 11.2. We define k (t)= £ h Ak h n (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. A n d we have a continuous utility function \{nx\ + x ). 2 • sup/j Ex\k (t)\ h < 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) = E | £ X + (tf,2) ) **(tf)} • \e- H^n,i) pt 2 2 A ( - ) 4 2 9 Since the chain is defined in the finite box G, we have Md,i) + (^, ) <Ci 2 2 2 . where C\ is a positive constant which is independent of h. This gives jf>-**Ai*(tf)J 4(0) < \c.Ex (4.30) As A i ' ( ^ ) —y 0, the sum in (4.30) is an approximation to the integral l poo 1= / Jo e-"*dt. The positivity of the p in our model guarantee that I < oo. Hence oo J2e- "At (^)<C pt h 2 n=0 where C is a positive constant independent of h. So 2 4(0) < \c c . x 2 Because C\ and C are independent of h, we have sup^ J^-(O) < oo. B y (3.23), 2 BViTpV (X) < sup J (0) h x h < oo. (4.31) h As everything in (3.22) is positive, by (4.31), we can say that sup Exi:^ oe-<\Ak \<oo h h su P / l = n E e" " IAfc£| < oo pt x < t n-|-l — =>• e sup^ E « =* c-"* s u p £ x | f c ( i ) l < _ p i x fc fc < because e "*» > e~ because E n A x , |A*£| = |* (*)|by (4.28) 0 0 _ pt h n+l — =>• sup^ Ex\k (t)\ h < oo because e~ pt So we have shown that sup^ Ex\k (t)\ h • i^n+i — '• } n I S a is a fixed number. < oo. 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 k to denote the optimal control for (3.23). In [16], Theorem 5.3 says h that k h converges weakly to some admissible control k(-) and the approximate Markov chain under the control k h 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.32) and (4.33) show that Jx(k(-)) — V(X). These observations say that our V (X) h the optimal k(X). This means that &(•) is an optimal control. and k (X) h (4.33) are approximations of the real V(X) and Chapter 5 Test Results In (2.1) and (2.5), we have parameters a i , a , 6n, 612, &215 ^22 c, fi and a. 2 5 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 dx (t) = (a + b2ix (t) 2 + b22X2(t))dt 1 2 + aidWtt) (5.1) + a dW(t) 2 with initial values (2.2) and cost function Jx(0) = E x | ^ ° ° 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, b n = - 1 , a = 0.1, b 2 2X = -1, b 22 = -3, a x = (0.03,0.01), a 2 = (-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, A n = -0.2, A 12 = -0.2, A 21 = 0.3 and A 22 = 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 v left line (-0.2, -0.2) horizontal middle line (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 = h \ — 0 . This means that the drift in the x\ direction is independent of 2 x and vice versa. Since the control is only exerted in the x\ direction, the decision whether to or 2 not to apply the control is not going to influence the process in the x direction. Consequently, 2 the control is not going to influence the cost arising from the x 2 direction. Hence we expect the region in which the control is active to be perpendicular to the x\ direction. exactly what we get from our numerical test; 6 2 2 = - 3 , [-0.2,0.3] ai = x ( 0 . 0 3 , 0 . 0 1 ) , tr [-0.2,0.3] with h 2 = (-0.01,0.02), = 0.005. We used a\ = 0 . 0 5 , &n = — 1 , a This is 2 ji = 0.1, a = 0.003 and p The control is shown in Figure 5.5. = 1 = 0.1, on the region 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 J (0) and J (0) x -0.2 x -0.15 -O.I -O.OS O : X = (x -0.2) lt COS 0.1 0.16 0.2 0.26 0.3 J (0) x -0,2 -0.1S and J (0) x -0.1 h = 0.01 J (0) and J (0) x -0.2 x -0.15 -0.1 -0.05 O 0.1 0.15 0.2 0.25 0.3 J (0) and J (0) x -0.2 x -0.15 -0.1 h = 0.01 J (0) x -0.2 -0.15 and J (0) x -0.1 —O.OS O 0.1 h = 0.01 O.OS 0.1 0.15 0.2 0.2S 0.3 -0.05 O : X = (a:i,0.05) 0.05 0,1 0.15 0.2 0.25 0.3 hr'= 0.005 J (0) : X = («i,0.3) 0.05 O -0.2) h = 0.005 : X = (si,0.05) 0.05 -0.05 :X = 0.15 0,2 0.25 x 0.3 2 -0.1 S and J (0) x -0.1 -0.05 O :X = 0.0S 0.1 0.1S h = 0.005 Figure 5.3: Comparison of different grid sizes on horizontal lines (x 0.3) u 0.2 0.25 0.3 Chapter 5. Test Results J (0) and j £ ( 0 ) : X = (-0.2, x ) x -0.2 < M ° ) and j £ ( 0 ) : X = (-0.2, z ) 2 -0.15 —O.I -0.05 O O.OS 0.1 0.15 0.2 2 0.3 0.25 .2 -0.16 -0.1 h = 0.01 J (0) -0.2 J (0) 2 -0.15 -0,1 -O.OB O O.OS 0.1 0.15 0.2 0.3 0.25 and J (0) x -0,2 x -0.15 -0.1 h = 0.01 J (0) x -O.Z -0.15 and J (0) x -0.1 -0.05 0 h = 0.01 -0.05 0 Jx(0) and J (0) 2 0.1 0.05 0.1 0.15 0.2 0,3 0.25 : X = (0.05, x ) 2 0.05 0.1 0.15 0,2 0.3 0.25 h = 0.005 : X = (0.3, x ) 0.05 O h = 0.005 and j £ ( 0 ) : X = (0.05, x ) x -O.OS 0.15 02 0.25 x O. -0.2 -0.15 -0.1 -0.05 0 : X = (0.3,x ) 2 0.05 0,1 0.15 h = 0.005 Figure 5.4: Comparison of different grid sizes on vertical lines 0,2 0.25 O. Chapter 5. Test Results 38 0.31 1 1 1 1 ^^^^5l»il^^^^^^^ o.25 iiiiiiiiiii 0.2 IIIIIIIIIII Jill 0.15 IIIIIIIIIII |§§I fijjjjj 0.1 Ijllllllll 0.05 | o IIIIIIIIIII BB -0.05 -0.1 IIIIIIIIIII IB HHHHHHHI -0.15 IIIIIIIIIII _Q2^ ^0.2 1 -0.1 -0.15 IK I -0.05 I 0 ' 0.05 X1 ' 0.1 0.15 0.2 0.3 0.25 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, a = 0.05 , 6 2 2i = —1, 622 = — 4, o\ = (0.03,0.01) and a = (-0.01,0.02) with p = 1, fj, = 0.1 and a = 0.002. We calculate the control 2 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 differential 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. Chapter 5. Test Results 39 0.25 0.2 0.15 0.1 -0.2 -0.15 -0.1 -0.05 0 0.05 XI 0.1 0.15 0.2 0.25 0.3 0.1 0.15 0.2 0.25 0.3 0.1 0.15 0.2 0.25 0.3 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.2i!iii!^ 0.15 1 -^0.2 -0.15 -0.1 -0.05 0 0.05 xi h = 0.005 -0.2 -0.15 -0.1 -0.05 O 0.05 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 -OAS -0.1 -0.05 0 40 0.1 0.05 h = 0.01 kh= 0.2 0.15 0.25 0.3 .2 -0.15 0.001 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 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 case 1 case 2 a\ 0.05 0.05 6n -1 -1 a-2 0.23 -0.09 &12 0.5 0.5 &21 &22 -2.5 2.5 -4 -3 o-u 0.03 0.03 012 Cr l 0.01 0.01 -0.01 -0.01 2 0~22 P 0.02 1 0.02 1 t 0.1 0.1 1 a 0.002 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\ dx^t) dt dx (t) dt 0.05 -x (t) 1 + 0.5x (t) 2 (5.3) 2 and ( 0 . 2 3 - 2.5 Xi(t) - 4x (t) 2 dx {t) 0.05 - z i ( i ) + 0.5 x (t) dt (5.4) dx (t) -0.09 + 2.5 x^t) - 3x (t). dt (5.3) and (5.4) are autonomous systems where the independent variable t does not appear x 2 2 2 explicitly. B y 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 ^0.2 -0.15 -0.1 -0.05 0 0.05 xl 41 0.1 0.15 0.2 0.25 0.3 -0 2 - 0 15 -0 1 case 1 - 0 05 0 0 05 xl 0 1 0.15 0.2 0.25 0.3 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. A n d 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\ , a , 621, b , c> o"i and o2 to give the model 2 2 22 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 Bill 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 x from the Canadian Consumer Price 2 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. Jan83 10 Jon94 20 Janes 30 40 JenBB SO JanB7 SO Janes Figure 6.1: Canadian treasury bill rates 44 TO DecBS Chapter 6. Data Fitting 45 115 10 jQn(34 Jan83 20 Jan85 30 AO SO JanB7 JanBB 60 JanSS 70 Doc88 Figure 6.2: Consumer price index 10 Jon84 Jan83 20 Jan85 30 AO JonSO SO Jan87 SO Jan88 70 D0088 Figure 6.3: The Bank of Canada rates 6.2 Estimation 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 Azi(f) = (a + b x (t) 1 11 l + b x (t))At 12 2 + a AW(t) 1 +cAk(t) (6.1) Ax {t) 2 = where AW(t) = \ZA~t(ei,e )' 2 (a + b x (t) + b x (t))At 2 21 1 22 2 + a AW(t) 2 and the e,-'s {i = 1,2) are independent random variables with standard normal distribution, because W(t) is an 7£ -valued standard Brownian motion. Now 2 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 X (t) is the continous process c that satisfies (2.1), then we have (cf. Kloeden and Platen [14], p.343) £ s u p \X (t) - X(t)\ < I(VA1 c t Chapter 6. Data Fitting 46 10 Jan64 Jon03 20 30 JonOS Jan86 AO Jon87SO SO Jan88 70 D0088 Figure 6.4: The inflation rates where K is a constant. This shows that X{t) is a strong approximation of order 1/2. A n d for each smooth function g, we have (cf. Kloeden and Platen [14], p.473) sup \E{g(X (t)) - g(X(t))}\ c t < K At g where K is a constant depending on g. This shows that X(t) is a weak approximation of order g 1. We use the notation xi(i) for £i(2,), Ax (i) x and k(i) for k{ti). If CPI(i) for Axi(ti), x (i) for x (ti), Ax (i) 2 2 stands for the consumer price index at time which is the inflation rate, by x (i) = (CPI(i) 2 - CPI(i - l))/CPI(i thus derived are shown in Figure 6.4. We derive Ax-i(i), Ax (i) 2 x\(i + 1) — Xi(i), Ax (i) 2 for we calculate Ax (ti) 2 x (i), 2 - 1). The inflation rates and Ak(i) using Axi(i) = = x (i + 1) — x {i) and Ak{i) — k{i) — k(i — 1) respectively. We use a 2 2 2 different way of calculating Ak(i) because it is to be a control variable, hence adapted. If we let y (i) = Axj(i)(j = 1,2), / ( t ) = Ak(i), 3 CXJ = o,-At(j = 1,2), 0 kj = b At(k,j = kj 1,2) and j j = 0kjVA~i(k,j = 1,2), then we can write (6.1) as k + Pi2X2(i) yi(0 = « i + Pnxi(i) 2/2(0 = a + /32ia;i(0 + / 22a;2(0 + (72i,722)(£i,£2) - + cf(i) + (711,7i2)(£i, £2)' 3 , (6.2) x / 2 (6.2) is a linear statistics model. We make different assumptions on Pi and c and compare the 2 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 i n case 3, we Chapter 6. Data Fitting case case case case 47 1 2 3 4 012 0 = 0 0 = 0 be estimated be estimated fix fix 012 is to /?12 is to 1 2 1 2 c fix c = 1 c is to be estimated fix c = 1 c is to be estimated Table 6.1: Four cases move f(i) to the other side of the equation and let y[(i) = yi(i) — f(i), then we have 2/i(0 = y (i) = 2 i + 0 n * i ( O + (3 x (i) + ( 7 n , 7 i 2 ) ( £ i , £ 2 ) ' a 12 2 (6.3) a2 + 0 2 i * i ( O + 022*2(O + (72i,722)(£i,£2)'. Since £\ and s are mutually independent standard normal variables, and are independent 2 of Xi and x , 2 we have the following equations if we let £i = 6 = 2/2 - OL - 0 i * i 2 y[ — a\ — fl\\X\ — (i\ x 2 2 and fl x : 2 22 2 E^ =0 E[(Zi ~ E^)( - E )] Xl Xl E[({i - E^)(x - Ex )} = E(^x ) 2 2 2 =0 =0 E£ 2 E[(h ~ E^)(x! - E )] Xl E[({ ~ Et )(x 2 ^ i = £(6*i) = 0 2 2 = £(6*i) = 0 - Ex )} = E(Z x ) 2 2 = 7?i + 7 i 2 2 2 =0 2 £ f l 6 = 711721 + 712722 E£, = 721 + 7 2 2 2 If we use "average" for the expectation, then we have N N 1 X ) ( y i ( 0 - « i - 0 n * i ( O - Pi2x (i)) 2 (6.4) =o N [ ( » i ( 0 - "1 - 011*l(O - 012*2(0) *!(*)] = N 0 (6.5) i=l N jj E M(i) i=l 1 N 2 E(^(*) ~ 4= 1 )9llXl(*) " 012*2(0) ^2(0] = 0 -Cti- 0 1 2 ~ 021*1(0 - 022*2(0) = 0 (6.6) (6.7) Chapter 6. Data Fitting 48 1 N J £ [(2/2(0 - "2 - 0 2 1 * l ( O - 0 2 2 * ( O ) * i ( O ] = (6-8) 0 2 i=i 1 X ^ [(2/2(0 - £ «2 - feariC*) - 0 2 2 * ( O ) * 2 ( O ] = 0 (6.9) 2 i=l 1 " ]yT E ( 2 / i ( 0 - « i - / ? i i * i ( 0 - & 2 Z ( 0 ) 2 2 = 7 u + 712 (6.10) Jf i£=l [ ( » l ( 0 - ^ 1 - 0 1 1 * l ( O - 0 1 2 * 2 ( 0 ) 1 8=1 (2/2(0 - « 2 - 0 2 1 * l ( O - 0 2 2 * 2 ( 0 ) ] = 711721 + 712722 N 1 (6.11) N E(2/2(0 - N «2 - 021*1(0 - 022* (O) 2 2 = 721 + 72 22 2 (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 , 0 i and / J . (6.10)-(6.12) give 2 2 2 2 the elements in the covariance matrix of £1 and £2- We assume that 712 = 7 i , then we can find 2 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 £ , that 2 is, the information about how £1 and £ correlate with each other, rather than the information 2 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 i and 6 2 22 are negative while &12 positive. This Chapter 6. Data Fitting ai 611 bi2 c &2 hi h2 Oil 012 = 0- l 2 022 equilibrium point 49 case 1 0.171 -1.805 0 1 0.971 -4.151 -13.974 0.0199 0.0005 0.0830 (0.095,0.041) case 2 0.116 -1.215 0 0.354 0.971 -4.151 -13.974 0.0156 0.0004 0.0830 (0.095,0.041) case 3 0.149 -1.705 0.307 1 0.971 -4.151 -13.974 0.0199 0.0005 . 0.0830 (0.095,0.041) case 4 0.084 -1.058 0.446 0.341 0.971 -4.151 -13.974 0.0152 0.0004 0.0830 (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 + 6 ) ± \/(b 22 u + 6 ) - 4(fen6 2 ~ 2 buhl) 2 22 2 The part under the square root is always less than |6n + 6221 because 6n6 2 22 — 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 6 i in our estimation. It might be that the data 2 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 x 2 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 02 .6 " 1 I 0 0 5 case 1 s is 1 0.25 0.3 0.2 0.25 0.3 case 2 0 1 s i 0.2 S 005 005 0 - 0 05 -0 1 - 0 15 -0.2^ 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 interest rates 0.15 0.2 0.25 0.3 case 4: p = 0.1 0.3 j : 0.25 tllltl lllllliil i.i i i i l l ! interest rates case 4: p = 0.05 0.3 0.25 0.2 -0.15 jllJIf llf?' _ Q 2 Li -0.2 I -0.15 ^^^^^^^^^^^^M I -0.1 I -O.05 I O I i 0.05 0.1 interest rates \nifi 11 II T ftfflffftiifrfrfrffiiftftiiihiii 0.2 0.25 0.3 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 01 | 0.0! ~°-0.2 -0.15 -0.1 -O.OS 0 0.05 0.1 0.15 mterest rales 0.2 0.25 0.3 case 4: p = 0 on region [-0.2,0.3] X [-0.2,0.3] ^ mimmm ^ mmm ;; ^^M wm 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 15 -02 -015 -0.1 -0.05 0 0.15 0.2 0.25 0.3 -0.2' -0.2 -0.15 case 4: a = 0.002 -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. Chapter 7 Conclusion 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 parameters and hence different control results. A t 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 ities p (X,Y) (X,Y G Gh) and time step At (X), h with transition probabil- that approximates the uncontrolled process, h then by the same argument as in the first part of Section 4.1, the value function V (X) h at point X for this Markov chain satisfies the following dynamic programming equation: = e- V (X) h ^J2p (X,Y)V (Y) pAth h h + ^( ixl Y (A.l) + x )At (X) 2 f h 2 2 Now our objective is to find p (X, Y) ( X , Y e G ) and At (X). h h h 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 i and 2 that (3.5) implies , x) -^1 d B ( -~ B + oxi where Bf(X) ( x n * + ^ ) - y f f l ) _ h 1 y ( * ) - n * - « * ) h % = 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" ) + 22 + (Bt + Bi)h + (B+ + B )h + a 2 = / M l l l - 2 12 a a _ 2 l 12 + + 7^[|«u - | « i 2 - a h B 2 i2 + a l ]il + , 2) V x h hB^]V(x 1 58 x - h,x ) 2 ph ]V( ,x ) 2 Xl 2 Appendix A. Derivation of the Transition Probabilities for the Diffusion Case ~ \at + ^[\ 22 a ~ \ \ 2 + hB%]V(x ,X + h) a 2 1 2 + p-[f«22 - | a i 2 - \o-i2 + hB ]V(x ,x 2 + x x + Wl i2 ( i a If we use Q (X) h a h - h) 2 ( i - h,x - v 2 - h,x + h) 2 2 h) x 2 + h, x - h) + ^\a^ V(xi V x + \(px\ 1 + ^2 t / ? 2 i 2 ^ ( i + h'-> 2 + ) a 59 2 + x\). (A.2) as defined in (3.17), then the coefficient of V(xi,x ) becomes 2 ^ m ^ + ) o ^ ) ( = 1 + _^L). ( A If we divided both sides of (A.2) by (A.3) and note that when h is small, ph /Q (X) 2 h ,, is small, we can say that Also, for a real number r, we have r V( ,x ) Xl = = 2 e s* { QJ^X) [|«H - 7 + r~ = \r\. Then (A.2) is approximately equivalent to + V(X) \ l2 ~ \ \2 + a + Q^X)f2 n a a |aj" 2 _ h | i2 + a B l ](1 X + > fe5f]V(a:i - V l) h X /i,z ) 2 + QH^f)[2 22 - 2 1 2 - \ \ 2 + hB$]V(x ,X a 1 + ^ X ) ^ 2 a + Qf^x)k i2 ( i V + \(lixl V a 2 + h,x + h)+ Q 7 ^ X ) | a i " F ( x i -h,x - x 2 2 + Qiqx)k i2 ( ~L a 2 ~ § i 2 ~ ^ i 2 + k B ~ ] V ( a : i , a : - h) 2 2 a + h) a a + M 2 x + xl)Q$ e &. = rj - h) + g^\ai V(xi 2 2 - h) h,x + h)} 2 (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 probabilities. 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), p ((x ,x ),(xi h 1 h),(h, -h). So we have + h,x )) + h E X,n^ n H (h, h),(—h, -h),(-h, 2 p ((x ,x ),(x - h 2 1 2 1 h,x )) 2 \ ° / + p ((x ,x ),(x ,x 1 V h 1 2 p ((x ,x ), 1 \ 1 -h (xi + h, x + h)) + | h 2 2 (-h\ 2 1 2 | p ((x!,x ), ( h 2 - h,x - Xl h)) 2 h I h - h)) h -h / + + 2 p ((x ,x ),(x ,x + h)) + h P ((x ,x ), (x-i - h,x + h)) + h 1 2 i k ^ 2 p ((x ,x ), ( h 1 2 Xl + h,x - h)) 2 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^-^, A^][A^-^ n (Atf,i - * & , „ A t f ) i n A^]' ( A t f - E Ati )(Ae , a fl (A&i fl - ^, Ae^)(A^, n where A £ ^ ( i = 1,2) denotes the 2 2& Atf ) iB |2 Xtn tl - n 2 E A^ ) x<n <2 (Af£ - 2 & A t f ) i2 iB 2 i2 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 p (A^ h) h 1 = = p ((x ,x ),(x 1 j l -/0 = / ( A ^ 1 0) = 2 1 2 p\( x\,x ),(x + /(A^ + h,x )) h 2 + h,x 2 + h)) (xi + h,x 2 - h)) 1 + P (( X!,x ), = p \ xi,x ),(xi - + P ( xi,x ),(xi - h,x - h)) + P { x\,x ),(xi -h,x + h)) = p { xi,x ),(x-i,x + P{ h 2 2 h 2 h 2 2 h 2 2 2 h h,x )) + h)) 2 x\,x ),(xi,x 2 - h)), 2 so E(A^ - £|Xi) = E(Ati f - = 0 p (A^ A 2 2 (EA& y A = 0) + h p (Ati h 2 tl = h) + (-h) V ( A # h tl = i X -h) -{B (X)At\X)f x = a At (X) + {B^X^At^X^ h n - Bl(X)(At (X)) h . 2 o(At»(X)) Similarly, we get E(A& - 2 = E A& ) 2 Xtn a At (X) 2 + \B (X)\At\X)h h 22 - B (X)(At (X)) 2 2 h . 2 v y ' And E(A^ - E Ati )(Ae , A Xtn tl E(AU,ihA^ ) - 2 (h • 0)p (( , - E Ae , ) h n 2 Xtn n 2 (EA& )(EA& ) tl 2 x ), ( x + h, x )) + (-h • 0)p ((xi,x ), h h Xl 2 a 2 2 (xi - h, x )) 2 + (0 • h)p ((x , X ), ( S i , X + h)) + (0 • ( - h ) ) ^ ^ Z ) , 532 - h)) h 2 X +(h • h)p ((x ,x ), (xi + h,x h 1 2 h 2 -(B!(X)At\X))(B (X)At (X)) h 2 + h)) + ((-h) • (-h))p (( ,x ), (x - h , x - h)) h 2 +((-h) • h)p ((x ,x ),(xi 1 2 2 - h,x Xl 2 + h)) + (h • (-h))p ((x ,x ),(x x h 2 1 2 t 2 + h,x 2 h)) Appendix B. Verification of Consistency Conditions = a At (X) + = a At (X) + h 12 h 21 (-B (X)B (X)(At (X)) ) h 1 o(At (X)). h Here we use the fact that a\ = a \. 2 62 2 These equations have shown that (3.21) is true. 2 2 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 /•OO Jx(0) = / 1 -ExixKt^e-^dt. (Cl) Now the problem of looking for the cost function is turned into looking for Ex[x\{t)). By I t o ' s lemma, d(x\(t)) = 2x (t)(a 1 + bux^tydt 1 + a dt + 2x (t)a dW(t) xl 1 1 (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 (0) + [ (2a x (t) + 2buxl(t) + a ) d f + 2 f (t)dW(t). Jo Jo 2 1 1 n Xl (C.3) If we take expectation on both sides, then we have E x\{t) f\2a x (t) + 2b xl{t) + a )dt + 2E 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: x = E xl(0) x +E x 1 1 E u X n f' (t)dW{t) Jo Xl X =0 because the stochastic integral is a martingale (cf. [12],pp.137-139). A n d if we switch the expectation and integral in (C.4), then we can write (C.4) as E x\(t) x = E xl(0) x + f E (2a (t) Jo x lXl 63 + 2b xl(t) n + a )dt, u Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution 64 or dE (xl(t)) = (2a E (x (t)) x 1 x + 2b E (x (t)) + a )dt. 2 1 11 x (C.5) u Similarly, dE {x\(t)) = x (2a Ex(x (t)) 2 + 2b Ex(x (t)x (t)) 2 21 1 + 2 2b E (x (t))dt 2 22 x 2 (C.6) dE (xi(t)x (t)) x = 2 (a E (xi(t)) 2 1 +b E (x (t)) x 2 + \{a 2 21 + (6 + a E (x (t)) x x U + b^Exix^x^t)) + a ))dt 21 (C.7) 12 dE (xi(t)) = ( dE (x (t)) = (a + b2iEx( (t)) + b22Ex(x (t)))dt x x 2 + b E (x (t)))dt fll 11 x (C.8) 1 2 Xl (C.9) 2 (C.5) through (C.9) is a hnear system of differential equations for Ex(xi(t)), Ex(x\{t)), Ex(x (t)) 2 and Ex{x\(t)x (t)) 2 with initial condition 2 Ex(xi(0)) = X l E (x (Q)) 2 = x Ex(xlm = xl Ex(xKO)) = x\ x Ex(x (t)), £x(xi(0)z (0)) 2 2 = (CIO) XxX 2 Starting from (C.8), we can solve (C.8), (C.9), (C.5), (C.7), (C.6) one by one for E (x (t)), x E (xl(t)), 2 x Ex(x (t)) 1 E ( (t)x (t)), x Xl and E (x (t)). 2 2 x E (x {t)), x The result is = (x + ^ ) e ^ - ^ b 1 On "ii C12 JP / / , u 021C11 E {x {t)) = 7- e b x 2 On - t a -b \C\ +( 7 ( 2 2 b \C\\ 7 7 - + x) + 2 On - (x (t)) 2 = e*»* + {x\ + 'n C31 02 a u -J 2 2 —O22 2 C22 C21 Ex 2 2 O22 °22 a -b \C\ 2 a i ° 2on C32 1 2 + ^ C23 i on ) e »t 2b + " n "Jj*" * " 1 -2on C*33 7 x Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution E ( (t)x (t)) x Xl = J^- 2 e »< + e^ +- 5 " b b C42 C41 C43 022 „ O i l + t>22 ^ 2 e On - 0 2 Oil -(On + 0 2 2 2 ) where = 1(^21 D = c ^ C i i + 01C21 + O21C31 D = «lC"22 D14 = 021C32 12 13 + O i 2 ) - a 2 C i 2 + O1C23 + O21C33 and Ex(4(t)) = - r - ^ r - e "'+ ^ + 6 Oil — ^°22^ C^iS • C51 ' 26 D + - + 7 1 C52 1 (ElL 1 + 2 1 22 e - C53 + '^?3 _ ^ 26 2 - J. On 2611 —OJ. 2622 OJ, J.022 OJ. ° 2 22 ^25 2 x 26 < 6 n u-*22 ^ x 2 ) e u 4 22 + 2 22 ^ e pu ~ 6 2 2 (6n+i»22)< ^ D 2 1 + ^-2622, y C55 C56 where -D21 = Z>22 = «22 + 2 a C 3 + 2621C45 2 2 a C i + 26 iC i , 2 2 2 D = 2a C D24 = D25 = , 2621C44. 23 2 2 2 2 + 26iC 4 42 2621C43 Note that we solve the differential equations when b\\ 7^ 622 and 6 n ^ 2 6 case for our test problem. 22 because this is Appendix C. Cost Function of the Uncontrolled Problem: the Exact Solution W i t h the solution of Ex{x {t)), 2 66 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( 2(t))e~ pt X is exponentially decaying. We have j /QN _ r 5 1 ' 2 -b +p l X K g K 11 C$2 -b 2 + P 2 C53 -2hi + P C54 -26 + P 22 C55 -bu - 6 22 + p Cjtf p' Bibliography [1] Akian, M . , Menaldi, J . L . , and Sulem, A . (1996). On an Investment-Consumption with Transaction Costs, working paper. Model [2] Bancora-Imbert, M . C . , Chow, P.L., and Menaldi, J . L . (1988). On the Numerical Approximation 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 Comparison 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: Bank Problem. Working paper. a Central [6] Fitzpatrick, B . G . , and Fleming, W . H . (1991). Numerical Methods for an Investment-Consumption Model. Math. Oper. Res. 16, No.4, 823-841. Oprimal [7] Fleming, W . H . , Zariphopoulou, T . (1991). An Optimal Investment/Consumption with Borrowing. Math. Oper. Res. 16, No.4, 802-822. Model [8] Hansen, L . P . (1982). Large Sample Properties of Generalized Method of Moments Estimators. 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 Equations. Springer, New York. 67 Bibliography 68 [15] Kushner, H . J . , and Dupuis, P . G . (1992). Numercial Methods for Stochastic Control Problems in Continuous Time. Springer, New York. [16] Kushner, H . J . , and Martins, L . F . (1991). Numerical Methods for Stochastic Singular Control Problems. S I A M J . Control and Optimization 29, No.6, 1443-1475. [17] Marsh, T . A . , and Rosenfeld, E . R . (1983). Stochastic Processes for Interest Rates and Equilibrium Bond Prices. J . Finance X X X V I 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 and Economic Policy edited by Fisher, S., U . of Chicago Press. Expectations [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 Control and Optimization 32, N o . l , 59-85. Models with Constraints. SIAM J.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical solution of a control problem in economics
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical solution of a control problem in economics Han, Bing 1996
pdf
Page Metadata
Item Metadata
Title | Numerical solution of a control problem in economics |
Creator |
Han, Bing |
Date Issued | 1996 |
Description | 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. |
Extent | 3787765 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-02-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079781 |
URI | http://hdl.handle.net/2429/4507 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
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
Cite
Citation Scheme:
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>
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