Construction and simulation of continuous time portfolio with given properties Sergey G. Shorokhov1 1 Peoples’ Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation Abstract A continuous-time version of portfolio management problem for a market with multiple stocks is under study. Asset allocation policy is determined implicitly as a solution to the system of ordinary differential equations for asset quantities. Right-hand sides of equations are derived explicitly in closed form from given portfolio properties using techniques of construction of differential equations by given integral manifold. Related issues of simulation for portfolio with desired properties are also addressed. We present the results of computer experiment for verification of specified portfolio property using constructed portfolio model in the form of the system of stochastic differential equations. Keywords Markov decision process, portfolio policy, stochastic process, SDE simulation, OpenCL framework 1. Introduction Modeling with stochastic differential equations is used extensively in many areas of modern research [1]. Stochastic differential equations (SDE) arise in modeling and simulation of various random dynamical systems in physical [2], chemical [3], biological [4], financial [4] and other sciences. Stochastic models based on SDE are an integral part of quantitative analysis of wireless channels in communications [5]. Modeling of signals and interference in communications requires construction of SDE with specified properties, such as given marginal probability density function and autocovariance function. Applications of SDE in other areas usually also imply construction of SDE with desired characteristics. We study asset allocation problem for a portfolio with specified properties using continuous- time continuous-state Markov decision process, determined implicitly by a system of differential equations. Our goal is to derive right-hand sides of equations by known portfolio properties and verify the required portfolio property in a computer experiment using simulation techniques. Workshop on information technology and scientific computing in the framework of the XI International Conference Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems (ITTMM-2021), Moscow, Russian, April 19–23, 2021 Envelope-Open shorokhov-sg@rudn.ru (S. G. Shorokhov) GLOBE https://esystem.rudn.ru/users/3295 (S. G. Shorokhov) Orcid 0000-0001-6835-4110 (S. G. Shorokhov) Β© 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) 33 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 2. Implicit Markov decision process for portfolio management Markov decision process [6] is a widely used approach to decision making, including applications in finance. Consider a market in which 𝑛 no dividend stocks are traded continuously and transaction cost and consumption are not taken into account. The stock price processes 𝑠𝑖 (𝑑) , 𝑖 = 1, 𝑛 generally satisfy the system of SDE [7]: π‘š 𝑑𝑠𝑖 = πœ‡π‘– (𝑑, 𝑠) 𝑑𝑑 + βˆ‘ πœŽπ‘–π‘— (𝑑, 𝑠) π‘‘π‘Šπ‘— , 𝑖 = 1, 𝑛, π‘š β©½ 𝑛, (1) 𝑗=1 where πœ‡π‘– (𝑑, 𝑠) are drift terms, πœŽπ‘–π‘— (𝑑, 𝑠) are volatility terms, π‘Š1 , ..., π‘Šπ‘š are independent standard Wiener processes, 𝑠 = (𝑠1 , ..., 𝑠𝑛 ). The initial stock prices at 𝑑 = 𝑑0 are assumed to be known: (0) 𝑠𝑖 (𝑑0 ) = 𝑠𝑖 > 0, 𝑖 = 1, 𝑛. (2) Let π‘₯𝑖 (𝑑) be the quantity (in units) of the 𝑖-th stock in the portfolio at time 𝑑. We assume that the quantity π‘₯𝑖 (𝑑) can be positive (long position), negative (short position) or zero (no position) and can take fractional values. The value process of the stock portfolio 𝑃 (𝑑) at time 𝑑 is determined by the formula 𝑛 𝑃 (𝑑) = βˆ‘ π‘₯𝑖 (𝑑) 𝑠𝑖 (𝑑) . (3) 𝑖=1 The share (weight) 𝑀𝑖 (𝑑) of the 𝑖-th stock in the portfolio 𝑃 (𝑑) is equal to π‘₯𝑖 (𝑑) 𝑠𝑖 (𝑑) π‘₯ (𝑑) 𝑠𝑖 (𝑑) 𝑀𝑖 (𝑑) = = 𝑛𝑖 , 𝑖 = 1, 𝑛. (4) 𝑃 (𝑑) βˆ‘π‘—=1 π‘₯𝑗 (𝑑) 𝑠𝑗 (𝑑) 𝑛 The weights 𝑀𝑖 (𝑑) are not all independent, because of the norming condition βˆ‘π‘–=1 𝑀𝑖 (𝑑) = 1. The state of the stock portfolio can be determined either by prices s = (𝑠1 , ..., 𝑠𝑛 ) and quantities x = (π‘₯1 , ..., π‘₯𝑛 ), or by prices s = (𝑠1 , ..., 𝑠𝑛 ), weights w = (𝑀1 , ..., 𝑀𝑛 ) and total portfolio value 𝑃 (𝑑). We assume that the portfolio state is given by price-quantity pair (s, x). The portfolio state (s, x) varies because of changes in stock prices s due to market situation or changes in stock quantities x due to policy (actions) of portfolio manager. Possible actions for the 𝑖-th stock in the portfolio are to hold the stock, to buy or to sell some units of the stock. As a result of actions for the time period △𝑑 > 0 change in the quantity of units of the 𝑖-th stock in the portfolio is equal to β–³π‘₯𝑖 = π‘₯𝑖 (𝑑 + △𝑑) βˆ’ π‘₯𝑖 (𝑑). If β–³π‘₯𝑖 > 0, then β–³π‘₯𝑖 units of the stock are purchased, if β–³π‘₯𝑖 < 0, then |β–³π‘₯𝑖 | units of the stock are sold, if β–³π‘₯𝑖 = 0, then π‘₯𝑖 (𝑑) units of the stock are being held in the portfolio. Asset allocation policy x (𝑑) may be explicit with quantities of stocks being functions of time and/or stock prices, i.e. x = x (𝑑, s). We assume that asset allocation policy x (𝑑) is determined implicitly as a solution of the system of differential equations 𝑑π‘₯𝑖 = 𝑓𝑖 (𝑑, s, x) 𝑑𝑑, 𝑖 = 1, 𝑛 (5) 34 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 with initial conditions (initial stock quantities) (0) π‘₯𝑖 (𝑑0 ) = π‘₯𝑖 , 𝑖 = 1, 𝑛. (6) Actually, the system (5) should be treated as a system of SDE with zero diffusion coefficients, because right-hand sides of system (5) depend on stochastic stock price processes 𝑠𝑖 (𝑑), following (1) with initial conditions (2). When functions 𝑓𝑖 (𝑑, s, x) in (5) are fixed, the portfolio management policy is fully determined by equations (1) and (5) with initial conditions (2) and (6), the portfolio dynamics depends only on the initial state of the portfolio (s(0) , x(0) ) at time 𝑑0 and is independent of the past (𝑑 < 𝑑0 ). Thus, we determine the portfolio policy by choosing various functions 𝑓𝑖 (𝑑, s, x) on the right-hand side of equations (5) and receive Markov decision process for portfolio management. 3. Construction of portfolio with given properties Construction of stock portfolio with desired characteristics and evaluation of its risk and performance using simulation and optimization plays extremely important role in modern financial industry [8]. The desired properties of the portfolio can be a given rate of return, a given risk metric (variance of return, value at risk, expected shortfall), given structural (country, industry, rating) ratios, etc. As we will show below, portfolio management policy may be based on construction of ordinary differential equations (ODE) by a given integral manifold. Originally, construction of ODE from a given integral curve was proposed by Yerugin [9] and later extended to dynamical systems of general nature by Galiullin [10], Mukharlyamov [11] and many other authors [12]. The method is widely used in investigations of controlled dynamical systems and allows to construct equations of motion from given properties of trajectories taking into account additional requirements, such as stability of the given manifold or optimality in some sense. The problem of construction of equations in the class of Ito SDE by known properties of motion was investigated by Tleubergenov [13, 14, 15]. We assume that price dynamics of stocks in the portfolio is described by SDE (1) with initial conditions (2) and quantities of stocks follow Markov portfolio policy, implicitly determined by equations (5) with initial conditions (6). The functions on the right-hand sides of equations (1) and (5) are assumed to be continuous in time 𝑑 and Lipschitz in portfolio state variables s, x. We consider portfolio properties in the form of 𝑙 equalities πœ”π‘˜ (𝑑, s, x) = 0, π‘˜ = 1, 𝑙, 𝑙 < 𝑛 (7) πœ• where Jacobian 𝑙 Γ— 𝑛-matrix ( πœ•x ) is of maximum rank 𝑙. Basically, the stochasticity of stock prices driven by SDE (1) leads to violation of equalities (7), so our goal is to find a portfolio policy of the form (5), which ensures that equalities (7) are satisfied on ensemble average, i.e. 𝔼 [πœ”π‘˜ (𝑑, s (𝑑) , x (𝑑))] = 0, π‘˜ = 1, 𝑙, βˆ€π‘‘ β©Ύ 𝑑0 , (8) 35 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 where 𝔼 denotes expected value (ensemble average) [16] and for a stochastic process 𝑋 (𝑑) +∞ 𝔼 [𝑋 (𝑑)] = ∫ π‘₯𝑓𝑋 (π‘₯, 𝑑) 𝑑π‘₯, (9) βˆ’βˆž where 𝑓𝑋 (π‘₯, 𝑑) is the probability density function of 𝑋 (𝑑) at time 𝑑. This statement of the problem differs from common approach to construction of stochastic differential equations by given integral manifold [13], when properties (7) are supposed to be satisfied exactly. But in portfolio management problem this implies restrictions on price equations (1), which is not relevant for financial market models. Stochastic differentials of πœ”π‘˜ can be calculated using the multidimensional Ito’s formula [17] 𝑛 𝑛 𝑛 𝑛 π‘š πœ•πœ”π‘˜ πœ•πœ” πœ•πœ” 1 πœ• 2 πœ”π‘˜ π‘‘πœ”π‘˜ = ( + βˆ‘ π‘˜ πœ‡π‘– + βˆ‘ π‘˜ 𝑓𝑖 + βˆ‘ βˆ‘ βˆ‘ πœŽπ‘–β„Ž πœŽπ‘—β„Ž ) 𝑑𝑑+ πœ•π‘‘ 𝑖=1 πœ•π‘ π‘– 𝑖=1 πœ•π‘₯𝑖 2 𝑖=1 𝑗=1 β„Ž=1 πœ•π‘ π‘– πœ•π‘ π‘— 𝑛 π‘š πœ•πœ”π‘˜ + βˆ‘ βˆ‘ πœŽπ‘–β„Ž π‘‘π‘Šβ„Ž , π‘˜ = 1, 𝑙. (10) 𝑖=1 β„Ž=1 πœ•π‘₯𝑖 Arbitrary (unknown) functions 𝑓𝑖 in (5) can be chosen in such a way that the following equalities hold true 𝑛 𝑛 𝑛 𝑛 π‘š 𝑙 πœ•πœ”π‘˜ πœ•πœ” πœ•πœ” 1 πœ• 2 πœ”π‘˜ + βˆ‘ π‘˜ πœ‡π‘– + βˆ‘ π‘˜ 𝑓𝑖 + βˆ‘ βˆ‘ βˆ‘ πœŽπ‘–β„Ž πœŽπ‘—β„Ž = βˆ’ βˆ‘ πœ†π‘˜π‘ž πœ”π‘ž , π‘˜ = 1, 𝑙, (11) πœ•π‘‘ 𝑖=1 πœ•π‘ π‘– 𝑖=1 πœ•π‘₯𝑖 2 𝑖=1 𝑗=1 β„Ž=1 πœ•π‘ π‘– πœ•π‘ π‘— π‘ž=1 where πœ†π‘˜π‘ž (𝑑, s, x) are arbitrary functions. Equations (11) can be regarded as a system of 𝑙 linear algebraic equations for 𝑛 unknown functions 𝑓𝑖 : 𝑛 πœ•πœ”π‘˜ βˆ‘ 𝑓𝑖 = βˆ’πœ‘π‘˜ , π‘˜ = 1, 𝑙, (12) 𝑖=1 πœ•π‘₯𝑖 where 𝑙 𝑛 𝑛 𝑛 π‘š πœ•πœ”π‘˜ πœ•πœ” 1 πœ• 2 πœ”π‘˜ πœ‘π‘˜ = βˆ‘ πœ†π‘˜π‘ž πœ”π‘ž + + βˆ‘ π‘˜ πœ‡π‘– + βˆ‘ βˆ‘ βˆ‘ πœŽπ‘–β„Ž πœŽπ‘—β„Ž . π‘ž=1 πœ•π‘‘ 𝑖=1 πœ•π‘ π‘– 2 𝑖=1 𝑗=1 β„Ž=1 πœ•π‘ π‘– πœ•π‘ π‘— Applying the methodology of Moore–Penrose pseudoinverse matrices [18] to system (12), we obtain the solution of (12) in the following matrix form: βˆ’1 πœ• + πœ• + πœ• πœ• + πœ•π‘‡ πœ• πœ•π‘‡ f = βˆ’( ) + [𝕀𝑛 βˆ’ ( ) ],( ) = ( ) , (13) πœ•x πœ•x πœ•x πœ•x πœ•x πœ•x πœ•x where πœ™ = (πœ‘1 , ..., πœ‘π‘™ )𝑇 , 𝕀𝑛 is the identity 𝑛 Γ— 𝑛-matrix, = (Ξ¦1 , ..., Φ𝑛 ) is a vector of arbitrary πœ• + πœ• functions. The pseudoinverse matrix ( πœ•x ) exists, because the Jacobian matrix πœ•x has the highest rank 𝑙. Substituting the corresponding expression from (12) for πœ™, we obtain the solution of (12) in the following final form πœ• + πœ• πœ• 1 πœ•2 πœ• + πœ• f = βˆ’( ) [Ξ› πœ” + + πœ‡ + tr (πœŽπ‘‡ 2 𝜎)] + [𝕀𝑛 βˆ’ ( ) ], (14) πœ•x πœ•π‘‘ πœ•s 2 πœ•x πœ•x πœ•x 36 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 where Ξ› = (πœ†π‘˜π‘ž ), tr is matrix trace operation. We additionally assume that for any π‘˜, π‘ž = 1, 𝑙 the random variables πœ†π‘˜π‘ž and πœ”π‘ž are indepen- dent, then, averaging the equalities (10) along ensemble of trajectories under the policy (14), we receive that 𝑙 𝑑𝔼 [πœ”π‘˜ ] = βˆ’ βˆ‘ 𝔼 [πœ†π‘˜π‘ž ] 𝔼 [πœ”π‘ž ] , π‘˜ = 1, 𝑙. (15) 𝑑𝑑 π‘ž=1 The system of ODE (15) admits the trivial solution 𝔼 [πœ”π‘˜ ] = 0, π‘˜ = 1, 𝑙, which implies that the portfolio with stock dynamics (1) and portfolio policy (5) admits portfolio properties (8). Thus, we receive the following statement on Markov portfolio policies with given portfolio properties. Markov asset allocation policy (5) with right-hand sides given by (14) admits the specified portfolio properties (8), provided that Ξ› = (πœ†π‘˜π‘ž (𝑑, s, x)) is an arbitrary 𝑙 Γ— 𝑙-matrix such that the random variables πœ†π‘˜π‘ž (𝑑, s (𝑑) , x (𝑑)) and πœ”π‘ž (𝑑, s (𝑑) , x (𝑑)) are independent, (𝑑, s, x) is an arbitrary column vector. Portfolio asset allocation policy (5) with right-hand sides given by (14) can be applied to portfolios containing stocks, currencies and commodities. By choosing arbitrary functions πœ†π‘˜π‘ž and Φ𝑖 it is possible to ensure the stability of portfolio policy or to find an optimal policy for given reward function and discount factor. Asset portfolio can include cash account, be self-financing, and incorporate other additional features. The value of portfolio with cash account and 𝑛 risky assets is equal to 𝑛 𝑃 (𝑑) = π‘₯0 (𝑑) + βˆ‘ π‘₯𝑖 (𝑑) 𝑠𝑖 (𝑑) , 𝑖=1 where π‘₯0 (𝑑) denotes the balance of cash account at time 𝑑. The portfolio is self-financing [17], if no external inflow or outflow of cash and stocks takes place and the purchase of a stock must be financed by cash on cash account or by sale of some stocks from the portfolio. The portfolio with cash account is self-financing, if [17] 𝑛 𝑑𝑃 = βˆ‘ π‘₯𝑖 𝑑𝑠𝑖 + π‘Ÿπ‘₯0 𝑑𝑑, 𝑖=1 where π‘Ÿ > 0 is the fixed interest rate of cash account. This implies that the balance of cash account π‘₯0 satisfies the following equation [19] 𝑛 𝑑π‘₯0 = (π‘Ÿπ‘₯0 βˆ’ βˆ‘ 𝑠𝑖 𝑓𝑖 ) 𝑑𝑑, 𝑖=1 where functions 𝑓𝑖 from (5) determine the portfolio policy. 4. Modeling self-financing portfolio with given structure We consider a self-financing portfolio with cash account and two risky assets (stocks), driven by geometric Brownian motion, and assume that the portfolio management policy is given in 37 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 implicit form (5). The portfolio dynamics is determined by the following system of SDE: 𝑑𝑠 = πœ‡1 𝑠1 𝑑𝑑 + 𝜎1 𝑠1 π‘‘π‘Š1 , ⎧ 1 βŽͺ 𝑑𝑠2 = πœ‡2 𝑠2 𝑑𝑑 + 𝜎2 𝑠2 π‘‘π‘Š2 , βŽͺ 𝑑π‘₯ = (π‘Ÿ π‘₯0 βˆ’ 𝑠1 𝑓1 (𝑑, s, x) βˆ’ 𝑠2 𝑓2 (𝑑, s, x)) 𝑑𝑑, (16) ⎨ 0 βŽͺ𝑑π‘₯1 = 𝑓1 (𝑑, s, x) 𝑑𝑑, βŽͺ βŽ©π‘‘π‘₯2 = 𝑓2 (𝑑, s, x) 𝑑𝑑. Here 𝑠1 , 𝑠2 are the stock prices, π‘₯0 is the balance of cash account, π‘₯1 , π‘₯2 are the quantities of stocks in the portfolio, s = (𝑠1 , 𝑠2 ), x = (π‘₯0 , π‘₯1 , π‘₯2 ), πœ‡1 , πœ‡2 are the instantaneous rates of return, 𝜎1 , 𝜎2 are the volatilities of stocks, π‘Š1 , π‘Š2 are independent standard Wiener processes, π‘Ÿ is the interest rate of cash account. We assume that the desired portfolio property is the condition that 80% of the portfolio assets is invested in risky assets (stocks). The value of the portfolio under consideration is equal to 𝑃 = π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 , and weights of cash and stocks are equal to π‘₯0 π‘₯1 𝑠1 π‘₯2 𝑠2 𝑀0 = , 𝑀1 = , 𝑀2 = . π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 The desired portfolio property can be expressed as the following relation between the stock weights 𝑀1 and 𝑀2 : π‘₯1 𝑠1 + π‘₯2 𝑠2 𝑀1 + 𝑀 2 = = 0.8, π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 which can be transformed into the equality πœ” ≑ βˆ’4π‘₯0 + π‘₯1 𝑠1 + π‘₯2 𝑠2 = 0. (17) Our goal is to determine the unknown functions 𝑓1 , 𝑓2 in (16) to guarantee the portfolio property (17) on ensemble average. According to Ito’s lemma the stochastic differential of πœ” is equal to π‘‘πœ” = (βˆ’4π‘Ÿπ‘₯0 + πœ‡1 π‘₯1 𝑠1 + πœ‡2 π‘₯2 𝑠2 + 5𝑠1 𝑓1 + 5𝑠2 𝑓2 ) 𝑑𝑑 + 𝜎1 𝑠12 π‘‘π‘Š1 + 𝜎2 𝑠22 π‘‘π‘Š2 . (18) We set the unknown functions 𝑓1 , 𝑓2 so, that the drift term at 𝑑𝑑 is equal to βˆ’π›Όπœ”, 𝛼 ∈ ℝ: βˆ’ 4π‘Ÿπ‘₯0 + πœ‡1 π‘₯1 𝑠1 + πœ‡2 π‘₯2 𝑠2 + 5𝑠1 𝑓1 + 5𝑠2 𝑓2 = βˆ’π›Όπœ”. (19) Equality (19) is the linear algebraic equation for the determination of unknown functions 𝑓1 , 𝑓2 , and the general solution of equation (19) is: 𝑠 𝑠 𝑓1 = βˆ’ 2 1 2 πœ‘ + 𝑠 2 πœ“ , 𝑓 2 = βˆ’ 2 2 2 πœ‘ βˆ’ 𝑠 1 πœ“ , (20) 𝑠1 + 𝑠2 𝑠1 + 𝑠 2 where πœ‘ = 0.2π›Όπœ” βˆ’ 0.8π‘Ÿπ‘₯0 + 0.2πœ‡1 π‘₯1 𝑠1 + 0.2πœ‡2 π‘₯2 𝑠2 , 𝛼 ∈ ℝ, πœ“ is an arbitrary function of 𝑑, s, x. 38 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 Thus, SDE system (16) for modeling prices and quantities of the portfolio with given property (17) takes the following form with an arbitrary constant 𝛼 ∈ ℝ and an arbitrary function πœ“: 𝑑𝑠 = πœ‡1 𝑠1 𝑑𝑑 + 𝜎1 𝑠1 π‘‘π‘Š1 , ⎧ 1 βŽͺ 𝑑𝑠2 = πœ‡2 𝑠2 𝑑𝑑 + 𝜎2 𝑠2 π‘‘π‘Š2 , βŽͺ βŽͺ𝑑π‘₯0 = (π›Όπœ” + 0.2π‘Ÿπ‘₯0 + 0.2πœ‡1 π‘₯1 𝑠1 + 0.2πœ‡2 π‘₯2 𝑠2 ) 𝑑𝑑, βŽͺ 𝑠 (21) βŽ¨π‘‘π‘₯1 = (βˆ’ 1 (0.2π›Όπœ” βˆ’ 0.8π‘Ÿπ‘₯0 + 0.2πœ‡1 π‘₯1 𝑠1 + 0.2πœ‡2 π‘₯2 𝑠2 ) + 𝑠2 πœ“) 𝑑𝑑, βŽͺ 𝑠1 + 𝑠22 2 βŽͺ βŽͺ βŽͺ𝑑π‘₯2 = (βˆ’ 𝑠1 (0.2π›Όπœ” βˆ’ 0.8π‘Ÿπ‘₯0 + 0.2πœ‡1 π‘₯1 𝑠1 + 0.2πœ‡2 π‘₯2 𝑠2 ) βˆ’ 𝑠1 πœ“) 𝑑𝑑. ⎩ 𝑠12 + 𝑠22 5. Simulation of self-financing portfolio with given structure To verify that the portfolio model (21) satisfies the portfolio property (17) we simulate the trajectories of system (21) for the time segment [𝑑0 , 𝑇]. 𝑇 βˆ’π‘‘ So we divide the segment [𝑑0 , 𝑇] into 𝑁 equal parts of length β„Ž = 𝑁 0 and simulate a discretized version of SDE (21). There exists a number of discretization schemes available, such as Milstein scheme [20] and stochastic version of Runge-Kutta methods [21, 22], but the most intuitive, easy to implement and common is Euler-Maruyama scheme [23]. Euler-Maruyama discretization of (21) gives us the following equations: 𝑠 = 𝑠 + πœ‡1 𝑠1,π‘˜ β„Ž + 𝜎1 𝑠1,π‘˜ βˆšβ„Žπ‘1,π‘˜ , ⎧ 1,π‘˜+1 1,π‘˜ βŽͺ𝑠 βŽͺ 2,π‘˜+1 = 𝑠2,π‘˜ + πœ‡2 𝑠2,π‘˜ β„Ž + 𝜎2 𝑠2,π‘˜ βˆšβ„Žπ‘2,π‘˜ , βŽͺπ‘₯ βŽͺ 0,π‘˜+1 = π‘₯0,π‘˜ + (π›Όπœ”π‘˜ + 0.2π‘Ÿπ‘₯0,π‘˜ + 0.2πœ‡1 π‘₯1,π‘˜ 𝑠1,π‘˜ + 0.2πœ‡2 π‘₯2,π‘˜ 𝑠2,π‘˜ ) β„Ž, 0.2π›Όπœ”π‘˜ βˆ’ 0.8π‘Ÿπ‘₯0,π‘˜ + 0.2πœ‡1 π‘₯1,π‘˜ 𝑠1,π‘˜ + 0.2πœ‡2 π‘₯2,π‘˜ 𝑠2,π‘˜ (22) ⎨π‘₯1,π‘˜+1 = π‘₯1,π‘˜ βˆ’ 𝑠1,π‘˜ β„Ž + 𝑠2,π‘˜ πœ“ β„Ž, βŽͺ 2 + 𝑠2 𝑠1,π‘˜ βŽͺ 2,π‘˜ βŽͺ 0.2π›Όπœ”π‘˜ βˆ’ 0.8π‘Ÿπ‘₯0,π‘˜ + 0.2πœ‡1 π‘₯1,π‘˜ 𝑠1,π‘˜ + 0.2πœ‡2 π‘₯2,π‘˜ 𝑠2,π‘˜ βŽͺπ‘₯2,π‘˜+1 = π‘₯2,π‘˜ βˆ’ 𝑠2,π‘˜ β„Ž βˆ’ 𝑠1,π‘˜ πœ“ β„Ž, 2 + 𝑠2 𝑠1,π‘˜ ⎩ 2,π‘˜ where 𝑠𝑖,π‘˜ = 𝑠𝑖 (π‘‘π‘˜ ), π‘₯𝑖,π‘˜ = π‘₯𝑖 (π‘‘π‘˜ ), π‘‘π‘˜ = 𝑑0 + π‘˜β„Ž, πœ”π‘˜ = βˆ’4π‘₯0,π‘˜ + π‘₯1,π‘˜ 𝑠1,π‘˜ + π‘₯2,π‘˜ 𝑠2,π‘˜ , 𝑍𝑖,π‘˜ are i.i.d. random variables with standard normal distribution, i.e. 𝑍𝑖,π‘˜ ∼ 𝒩 (0, 1). Basically, simulation of SDE is a resource-intensive application, so to speed up the simulation one has to enable GPU acceleration. PyOpenCL [24] is a programming environment for access to OpenCL parallel computation framework [25] from Python language. To apply PyOpenCL environment for portfolio simulation, one has to follow the typical sequence of steps: 1. Import the OpenCL API, obtain an OpenCL platform and a GPU device id: import pyopencl as cl plat = cl.get_platforms()[0] dev = plat.get_devices()[2] 39 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 2. Initialize a context for the selected GPU device, create a Queue object (with profiling enabled to track computation time), create memory buffers on GPU for input and output data (here p , q and s_gpu are NumPy arrays): opencl_context = cl.Context(devices=[dev]) command_queue = cl.CommandQueue(opencl_context, properties=cl.command_queue_properties.PROFILING_ENABLE) p_buffer = cl.Buffer(opencl_context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=p) q_buffer = cl.Buffer(opencl_context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=q) s_buffer = cl.Buffer(opencl_context, cl.mem_flags.WRITE_ONLY, s_gpu.nbytes) 3. Store the source code of C functions for GPU execution in Python strings (sde_step_src contains source code of one step of SDE simulation according to (22), simulate_sde_src contains source code of main kernel function, being called from Python): sde_step_src = """ // One step of SDE simulation static void sde_step(__global double *p, double *x, double W) { double s1n, s2n, x0n, x1n, x2n; x[5] = -0.8*x[2] + 0.2*x[3]*x[0] + 0.2*x[4]*x[1]; // omega s1n = p[P_M1]*x[0]*p[P_DT] + p[P_S1]*x[0]*p[P_SQRDT]*W; s2n = p[P_M2]*x[1]*p[P_DT] + p[P_S2]*x[1]*p[P_SQRDT]*W; x0n = (p[P_A]*x[5] + 0.2*p[P_R]*x[2] + 0.2*p[P_M1]*x[3]*x[0] + 0.2*p[P_M2]*x[4]*x[1])*p[P_DT]; x1n = -x[0]*(p[P_A]*x[5] - 0.8*p[P_R]*x[2] + 0.2*p[P_M1]*x[3]*x[0] + 0.2*p[P_M2]*x[4]*x[1])*p[P_DT]/(x[0]*x[0]+x[1]*x[1]); x2n = -x[1]*(p[P_A]*x[5] - 0.8*p[P_R]*x[2] + 0.2*p[P_M1]*x[3]*x[0] + 0.2*p[P_M2]*x[4]*x[1])*p[P_DT]/(x[0]*x[0]+x[1]*x[1]); x[0] += s1n; x[1] += s2n; x[2] += x0n; x[3] += x1n; x[4] += x2n; }""" simulate_sde_src = """ // GPU kernel function to simulate trajectories of SDE __kernel void simulate_sde(__global double *p, __global int *q, __global double *s) { // Indexing the current element (trajectory) to process int i = get_global_id(0); // Pointer to the i'th row of s (output) __global double *s_i = &s[i*6]; // 6 values in output // Simultaneous calculation of SDE trajectories within OpenCL kernel double x[6]; // placeholders for s1,s2,x0,x1,x2,omega 40 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 int q_i = (int) q[i], *seed = &q_i; const double r4_pi = 3.141592653589793; double v1,v2; for (int j = 0; j < 5; j++) { x[j] = p[P_S10+j]; } for (int j = 0; j < (int) p[P_N]/2.; j++) { // 2 steps of SDE v1 = r8_uniform_01 ( seed ); v2 = r8_uniform_01 ( seed ); sde_step( p, x, sqrt(-2.0*log(v1))*cos(2.0*r4_pi*v2) ); sde_step( p, x, sqrt(-2.0*log(v1))*sin(2.0*r4_pi*v2) ); } for (int j = 0; j < 6; j++) { s_i[j] = x[j]; } }""" 4. Compile the GPU C program from source: opencl_program = cl.Program(opencl_context, sde_hdr_src+sde_rng_src+sde_step_src+simulate_sde_src).build() 5. Enqueue the SDE simulation program on the GPU device and wait until the program completion: event = opencl_program.simulate_sde(command_queue, s_gpu.shape, None, p_buffer, q_buffer, s_buffer) event.wait() 6. Get back the output data from GPU memory into NumPy array s_gpu : cl.enqueue_copy(command_queue, s_gpu, s_buffer).wait() Function r8_uniform_01 is an implementation of uniform random number generator [26], C++ versions of other random number generators may be found in [27]. Standard normally distributed random variables 𝑍𝑖,π‘˜ are implemented with Box–Muller transform [28]. Simulation of discretized equations (22) is performed using computer program in Python with the following parameters: 𝑠1,0 = 0.9, 𝑠2,0 = 1.1, πœ‡1 = 0.3, 𝜎1 = 0.4, πœ‡2 = βˆ’0.1, 𝜎2 = 0.3, 𝑑0 = 0., 𝑇 = 1. π‘₯0,0 = 90., π‘₯1,0 = 110., π‘₯2,0 = 100., π‘Ÿ = 0.05, 𝛼 = 1., πœ“ ≑ 0, 𝑁 = 100. The results of the computer experiment are shown in Fig. 1. Green line on the last plot with stock weights shows that the total share of stocks 𝑀1 + 𝑀2 tends to the value of 0.8, which is the target property of the portfolio under consideration. Thus, constructed equations (22) generate trajectories with empirical portfolio property close to the given property (17). 41 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 Figure 1: Simulation of portfolio with cash and two stocks and given structural ratio 𝑀1 + 𝑀2 = 0.8 Acknowledgments The research was funded by RFBR, grant No. 19-08-00261. References [1] E. Allen, Modeling with ItΓ΄ Stochastic Differential Equations, Mathematical Modelling: Theory and Applications, Springer Netherlands, 2007. doi:10.1007/978- 1- 4020- 5953- 7 . [2] N. V. Kampen, Stochastic differential equations, Physics Reports 24 (1976) 171–228. doi:10.1016/0370- 1573(76)90029- 6 . [3] R. P. King, Applications of stochastic differential equations to chemical-engineering problems: An introductory review, Chemical Engineering Communications 1 (1974) 221–237. doi:10.1080/00986447408960433 . [4] C. A. Braumann, Introduction to Stochastic Differential Equations with Applications to Modelling in Biology and Finance, Wiley, 2019. doi:10.1002/9781119166092 . [5] S. Primak, V. Kontorovich, V. Lyandres, Stochastic Methods and Their Applications to Communications, Wiley, 2004. doi:10.1002/0470021187 . 42 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 [6] R. Bellman, A markovian decision process, Journal of Mathematics and Mechanics 6 (1957) 679–684. doi:10.1512/iumj.1957.6.56038 . [7] X. Y. Zhou, G. G. Yin, Markowitz’s mean-variance portfolio selection with regime switching: A continuous-time model, SIAM J. Control. Optim. 42 (2003) 1466–1482. doi:10.1137/ S0363012902405583 . [8] D. A. Pachamanova, F. J. Fabozzi, Portfolio Construction and Analytics, John Wiley & Sons, Inc, 2016. doi:10.1002/9781118656747 . [9] N. Erugin, Construction to all set of differential systems having a given invariant curve, Prikladnaia matematica i mehanika 16 (1952) 659–670. [10] A. S. Galiullin, Certain problems in the design of programmed-motion systems, in: Symposia on Theoretical Physics and Mathematics 8, Springer US, 1968, pp. 185–192. doi:10.1007/978- 1- 4684- 7721- 4_16 . [11] R. G. Mukharlyamov, On the construction of differential equations of motion of con- strained mechanical systems, Differential Equations 39 (2003) 369–380. doi:10.1023/a: 1026021701825 . [12] R. G. Mukharlyamov, M. I. Tleubergenov, Control of system dynamics and constraints stabi- lization, in: Communications in Computer and Information Science, Springer International Publishing, 2017, pp. 431–442. doi:10.1007/978- 3- 319- 66836- 9_36 . [13] M. I. Tleubergenov, An inverse problem for stochastic differential systems, Differential Equations 37 (2001) 751–753. doi:10.1023/a:1019285119532 . [14] M. I. Tleubergenov, On the inverse stochastic reconstruction problem, Differential Equations 50 (2014) 274–278. doi:10.1134/s0012266114020165 . [15] M. Tleubergenov, and G.T. Ibraeva and, On inverse problem of closure of differential systems with degenerate diffusion, Eurasian Mathematical Journal 10 (2019) 93–102. doi:10.32523/2077- 9879- 2019- 10- 2- 93- 102 . [16] O. C. Ibe, Markov Processes for Stochastic Modeling, 2nd edition ed., Elsevier, 2013. doi:10.1016/c2012- 0- 06106- 6 . [17] T. BjΓΆrk, Arbitrage Theory in Continuous Time, 4th edition ed., Oxford University Press, 2019. doi:10.1093/oso/9780198851615.001.0001 . [18] G. Wang, Y. Wei, S. Qiao, Generalized Inverses: Theory and Computations, Springer Singapore, 2018. doi:10.1007/978- 981- 13- 0146- 9 . [19] S. Shorokhov, Management of Financial Asset Portfolios, RUDN University, 2013. [20] G. N. Mil’shtejn, Approximate integration of stochastic differential equations, Theory of Probability & Its Applications 19 (1975) 557–562. doi:10.1137/1119062 . [21] P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer Berlin Heidelberg, 1992. doi:10.1007/978- 3- 662- 12616- 5 . [22] M. N. Gevorkyan, T. R. Velieva, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, Stochas- tic Runge–Kutta software package for stochastic differential equations, in: Dependability Engineering and Complex Systems, Springer International Publishing, 2016, pp. 169–179. doi:10.1007/978- 3- 319- 39639- 2_15 . [23] G. Maruyama, Continuous markov processes and stochastic equations, Rendiconti del Circolo Matematico di Palermo 4 (1955) 48–90. doi:10.1007/bf02846028 . [24] A. KlΓΆckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, A. Fasih, PyCUDA and PyOpenCL: A scripting-based approach to GPU run-time code generation, Parallel Computing 38 (2012) 43 Sergey G. Shorokhov CEUR Workshop Proceedings 33–44 157–174. URL: https://doi.org/10.1016%2Fj.parco.2011.09.001. doi:10.1016/j.parco.2011. 09.001 . [25] B. Gaster, L. Howes, D. R. Kaeli, P. Mistry, D. Schaa, Heterogeneous Computing with OpenCL, Elsevier, 2012. doi:10.1016/c2011- 0- 69669- 3 . [26] J. Burkardt, UNIFORM – a uniform random number generator (RNG), 2019. URL: https: //people.sc.fsu.edu/~jburkardt/c_src/uniform/uniform.html. [27] M. N. Gevorkyan, A. V. Demidova, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastianov, I. M. Gostev, Pseudo-random number generator based on neural network, in: CEUR Workshop Proceedings, volume 2267, 2018, pp. 568–572. [28] G. E. P. Box, M. E. Muller, A note on the generation of random normal deviates, The Annals of Mathematical Statistics 29 (1958) 610–611. doi:10.1214/aoms/1177706645 . 44