Solving High-Dimensional Dynamic Portfolio Choice Models with Hierarchical B-Splines on Sparse Grids

Schober, Peter; Valentin, Julian; Pflüger, Dirk Article — Published Version Solving High-Dimensional Dynamic Portfolio Choice Models with Hierarchical B-Splines on Sparse Grids Computational Economics Provided in Cooperation with: Springer Nature Suggested Citation: Schober, Peter; Valentin, Julian; Pflüger, Dirk (2021) : Solving High-Dimensional Dynamic Portfolio Choice Models with Hierarchical B-Splines on Sparse Grids, Computational Economics, ISSN 1572-9974, Springer US, New York, NY, Vol. 59, Iss. 1, pp. 185-224, https://doi.org/10.1007/s10614-020-10061-x This Version is available at: https://hdl.handle.net/10419/287127 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Vol.:(0123456789) Computational Economics (2022) 59:185–224 https://doi.org/10.1007/s10614-020-10061-x 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models withHierarchical B‑Splines onSparse Grids PeterSchober1 · JulianValentin2· DirkPflüger2 Accepted: 15 October 2020 / Published online: 4 January 2021 © The Author(s) 2020 Abstract Discrete time dynamic programming to solve dynamic portfolio choice models has three immanent issues: firstly, the curse of dimensionality prohibits more than a handful of continuous states. Secondly, in higher dimensions, even regular sparse grid discretizations need too many grid points for sufficiently accurate approximations of the value function. Thirdly, the models usually require continuous control variables, and hence gradient-based optimization with smooth approximations of the value function is necessary to obtain accurate solutions to the optimization problem. For the first time, we enable accurate and fast numerical solutions with gradient-based optimization while still allowing for spatial adaptivity using hierarchical B-splines on sparse grids. When compared to the standard linear bases on sparse grids or finite difference approximations of the gradient, our approach saves an order of magnitude in total computational complexity for a representative dynamic portfolio choice model with varying state space dimensionality, stochastic sample space, and choice variables. Keywords Curse of dimensionality· Dynamic portfolio choice· Discrete time dynamic programming· Gradient-based optimization· Spatially adaptive sparse grids· Hierarchical B-splines * Peter Schober [email protected] Julian Valentin [email protected] Dirk Pflüger dirk.pflueg[email protected]art.de 1 Finance Department, Goethe University Frankfurt, Theodor-W.-Adorno-Platz 3, 60323FrankfurtamMain, Germany 2 Institute forParallel andDistributed Systems, University ofStuttgart, Universitätsstraße 38, 70569Stuttgart, Germany 186 P.Schober et al. 1 3 1 Introduction A common approach to solve dynamic portfolio choice models in discrete time is dynamic programming, iterating over the value function backwards in time. Starting from the known value function at final time T, the value function is approximated on a state space grid, assuming that the state space is continuous. To determine the next iterate of the value function at each grid point at time T−1 we have to solve an optimization problem that depends on the value function of the previous iterate at time T. When a tensor product approximation is used, this approach suffers from the curse of dimensionality as the number of grid points of the approximation grows exponentially with the dimensionality of the state space. In addition, solving for the current value function iterate at a grid point relies on an accurate solution of the underlying optimization problem. When the portfolio choice is continuous, e.g., choosing the investment amount in stocks, bonds, etc., the computation of the optimal solution can be greatly accelerated by gradient-based optimization routines if the gradient of the objective function is available. Recently, sparse grids have been successfully employed to break the curse of dimensionality in high-dimensional dynamic models (Brumm and Scheidegger 2017; Judd et al. 2014; Schober 2018; Winschel and Krätzig 2010).1 A standard d-dimensional tensor product grid on the unit hypercube [0, 1]d with mesh size 2−n , n∈ℕ , and no points on the boundary contains 2n −1 grid points per coordinate direction and thus O(2nd) points in total, growing exponentially with the dimensionality d. In contrast, a regular sparse grid with the same mesh size contains only O(2nnd−1) points. The error of the sparse grid approximation of a function with homogeneous boundary conditions using piecewise linear basis functions is O(2−2nnd−1) with respect to the L2 and L∞ norm if the approximated function has bounded mixed second derivatives (Bungartz and Griebel 2004; Zenger 1991). This is only slightly worse than the corresponding error O(2−2n) for the case of full tensor product grids. In higher dimensions, even regular sparse grids need too many grid points for a sufficiently accurate approximation when solving high-dimensional dynamic models (Brumm and Scheidegger 2017). Fortunately, for approximations in the standard piecewise linear basis, the hierarchical structure of sparse grids allows for spatially adaptive refinement of the grid by inserting the 2d children of only certain leaves in the hierarchical structure. Spatially adaptive refinement has successfully been employed to solve high-dimensional dynamic models by Brumm and Scheidegger (2017) and Schober (2018). Unfortunately, approximations of the value function using the standard piecewise linear basis are not continuously differentiable and, hence, have discontinuous gradients. This poses a problem to gradient-based optimization techniques, which rely on a twice continuously differentiable approximation of the objective function to ensure convergence (Schober 2018). 1 See Cai (2019) for a discussion on alternative approximation methods for dynamic economic models that also do not suffer from the curse of dimensionality. 187 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… Global polynomial approximations have shown to work well with value function iteration and continuous choices for solving dynamic economic models (Cai and Judd 2015; Judd etal. 2014) as they are globally smooth. Smolyak’s formula can be used to construct sparse grid approximations (Barthelmann etal. 2000) on global polynomial bases, which can be refined adaptively with regard to specific dimensions of the state space (Judd etal. 2014) and with regard to the hierarchical surpluses, i.e., locally adaptively (Stoyanov 2017). Value function iteration with the use of gradient information to approximate the value function more accurately with global polynomials is also possible (Cai and Judd 2015). However, B-splines are much more flexible than global polynomials (Valentin and Pflüger 2016; Valentin 2019). While global polynomial approximations are bound to certain grid structures to avoid Runge’s phenomenon or similar issues, B-spline basis functions can be employed on any nested spatially adaptive grid hierarchy. They allow for simultaneous localand degree-adaptive refinement (hpadaptivity), implying that one could use a smaller or larger mesh size and/or degree of the B-spline basis functions in certain regions of the state space, e.g., to resolve kinks. In addition, the local basis functions are faster to evaluate than the conventional global polynomial basis functions. Approximations with B-splines of cubic degree (or higher) are twice continuously differentiable, and readily supply smooth and explicit approximations of both, the value function and the gradient. Compared to approximating the derivatives with finite differences, the optimization is not only more accurate but also significantly faster, especially when the number of optimization variables is large (Valentin 2019). B-splines have thus proven useful for computing numerical solutions to numerous dynamic models when finding the root of the gradient is required (Chu etal. 2013; Habermann and Kindermann 2007; Judd and Solnick 1994; Philbrick and Kitanidis 2001). In total, three issues with discrete time dynamic programming for dynamic portfolio choice models with continuous choices emerge: the curse of dimensionality, the lack of spatial adaptivity, and the lack of continuous gradients. It is apparent that current economic literature deals with these issues only in isolation, e.g., by combining sparse grids with global polynomial basis functions, or using sparse grids with non-smooth local linear basis functions to allow for spatial adaptivity. These approaches are hence computationally inefficient in accurately solving high-dimensional dynamic portfolio choice models or any high-dimensional dynamic economic model that requires smooth approximations or gradient-based optimization. This paper is the first to address all of these issues at once by combining hierarchical B-splines with sparse grids to approximate the value function and its gradient. Thus, we enable accurate and fast numerical solutions using gradient-based optimization while still allowing for spatial adaptivity (Pflüger 2010; Valentin and Pflüger 2016). The hierarchical grid structure allows us to develop an algorithm that uses the local adaptivity similar to Brumm and Scheidegger (2017) and Schober (2018), but interpolates the value function and its gradient with a B-spline basis. Therefore, we create a sparse grid for the value function, for which we interpolate the value function in the piecewise linear basis. We then refine the grid using the standard hierarchical surplus-volume-based refinement 188 P.Schober et al. 1 3 criterion. Finally, we interpolate the value function with hierarchical B-spline basis functions on the spatially adaptively refined sparse grid. We focus our study on the numerical accuracy of our approach. Therefore, we choose a dynamic portfolio choice model with multiple stocks, one bond, and consumption. For buying and selling the stocks, linear transaction costs have to be deducted. The resulting optimization problem is high-dimensional in terms of the state space, stochastic sample space, and choice variables. Hence, this problem is especially suited for a complexity analysis. At the same time, this model is similar to models from a vast strand of state-of-the-art literature on dynamic portfolio choice, e.g., Barberis and Huang (2009), Cocco etal. (2005), De Giorgi and Legg (2012), Horneff etal. (2010), Horneff etal. (2008), Hubener etal. (2016), Hubener etal. (2014), Inkmann etal. (2011). Consequently, our approach can be generalized to a broad class of dynamic portfolio choice models with only minor modifications. Dynamic portfolio choice models with transaction costs have been studied economically, e.g., by Abrams and Karmarkar (1980), Kamin (1975), Liu and Loewenstein (2002), Magill and Constantinides (1976), and extensively numerically by Cai (2009), Cai and Judd (2010), Cai etal. (2015), Cai etal. (2020). The latter report computational times and economic solutions for higher-dimensional transaction costs problems. They employ polynomial interpolation with only few polynomial nodes and parallelization to solve these problems with and without consumption using value function iteration in discrete time. Cai etal. (2020) present convergence results and computational times for the three-dimensional transaction costs problem with consumption and numerical errors for the four-dimensional problem using complete Chebyshev polynomials to approximate the value function. However, only we have employed spatially adaptive sparse grids to the transaction costs problem, which induces optimization on continuous choices (Schober 2018). We also apply local adaptivity to compute the optimal policies from the solution to the underlying optimization problem as earlier suggested by us (Schober 2018) and Brumm and Grill (2014). A complexity analysis reveals that cubic B-splines save more than one order of magnitude in computational effort compared to the state of the art with the linear basis (Brumm and Scheidegger 2017) and/or finite difference approximations of the gradient on regular sparse grids. Using spatially adaptive refinement of the optimal policy, we solve the problem for up to five dimensions where highly accurate solutions on regular sparse grids would require hundreds of thousands of grid points, and are thus no longer suitable. Here, spatially adaptive refinement allows a comparably low base resolution in the solution process of the value function, and, in a second step, adds grid points in the optimal policies where required. By this, we obtain low reported unit-free Euler equation errors for the transaction costs problem. The rest of this article is structured as follows: In Sect.2, we define the general class of dynamic portfolio choice models for which our approach is applicable. Section3 introduces hierarchical B-splines on spatially adaptive sparse grids, leading to the definition of hierarchical weakly fundamental not-a-knot splines. Algorithms for solving dynamic portfolio choice models with B-splines on spatially adaptive sparse grids are discussed in Sect.4. We analyze the complexity and demonstrate the numerical 189 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… accuracy of our approach solving the transaction costs problem in Sect.5 before concluding in Sect.6. 2 Discrete Time Dynamic Portfolio Choice Models We consider discrete time dynamic portfolio choice models with finite time horizon T in which the investor seeks to maximize additive expected life-time utility u from consumption ct : Here, she has mp continuous choices pt∈𝛹⊂ℝmp (e.g., investment amounts in stocks and bonds) with respect to d continuous states xt∈𝛺⊂ ℝ d (e.g., current financial wealth or labor income) she can reside in at time t. In addition, the transition from state xt to xt+1 does not only depend on her choices and her state, but may also be subject to m𝜻 random shocks 𝜻t∈Z (such as stock returns or labor income shocks), which are drawn from the sample space Z⊂ℝm𝜻 . The random variable ft∶𝛹×𝛺×Z → 𝛺 , (pt,xt,𝜻t) ↦ xt+1 , then describes the continuous state dynamics between t and t+1 . We denote by 𝜌<1 the subjective time discount factor, and we assume the utility function to be of Constant Relative Risk Aversion type with risk aversion 𝛾>1 : It is also straightforward to include discrete choices (compare also, e.g., Brumm and Scheidegger 2017) in the trivial way and discrete states (Schober 2018) in this model setup. Furthermore, the model can be generalized to further utility functions, e.g., to Epstein and Zin (1989) utility and to utility functions with a narrow framing component (Barberis and Huang 2009). In this paper, we disregard these modeling choices purely for simplicity. By the Bellman principle(Bellman 1954), this utility maximization problem can be reformulated in terms of the value function jt for t=0, …,T with known terminal utility v subject to the mg possibly non-linear inequality constraints gt∶𝛹×𝛺 →ℝ mg : For two vectors a,b , we define a≥b if ai≥bi for all i (and “ ≤ ” analogously). The corresponding expected value of jt+1 is (1) 𝔼 0 [T ∑ t=0 𝜌tu(ct(pt,xt)) ]. (2) u (ct)= 1 1−𝛾 c1−𝛾 t . (3a) jt(xt)=max p t { u(ct(pt,xt)) + 𝜌𝔼t [ jt+1(ft(pt,xt,𝜻t)) ]} ,t<T , (3b) jT(xT)=v(xT), (3c) gt(pt,xt)≥0. 190 P.Schober et al. 1 3 Here, 𝛷t(⋅|xt) denotes the conditional distribution of 𝜻t . Note that we only need the value function at time t+1 to determine the value function at time t via maximization. To numerically solve the optimization problem(3), discrete-time dynamic programming iterating over the value function is common (Judd 1998; Rust 2008). Therefore, the value function jt is restricted to a finite grid on the (truncated) state space with Nt grid points x(k) t , k=1, …,Nt . The value function values in between grid points are interpolated by with basis functions 𝜑k and coefficients 𝛼k , which are chosen in such a way that the interpolant j S t fits the known function values at all grid points x(k) t . Beginning with the known final solution at time T, the Bellman equation is solved backwards in time until the value function is computed for each grid point at each t=T−1, …,0 . For t<T , let us define the interpolant of the objective function of the maximization in Eq. (3a) for grid point x(k) t by: The maximization of this target function(6) with respect to pt can then be performed using Sequential quadratic programming (SQP) routines, see "AppendixA.1". To compute the expectation with respect to d𝛷 , numerical integration can be used if the conditional distributions 𝛷 (⋅ | x (k) t) are known.2 3 Hierarchical B‑Splines onSparse Grids As discussed in Sect. 1, hierarchical B-splines on sparse grids provide numerous advantages over other basis choices, especially in the context of optimization (Valentin 2019; Valentin and Pflüger 2016). In addition, they allow for spatially adaptive refinement by applying the standard surplus-based refinement criterion. To compute the coefficients of the B-spline approximation, usually a computationally expensive linear system has to be solved. Therefore, we determine the underlying grid structure by applying the surplus-based refinement criterion on the piecewise linear basis and interpolating with B-splines on the resulting grid. As proven in our previous work (4) 𝔼 t [ jt+1(ft(pt,xt,𝜻t)) ] ∶= ∫ Z jt+1(ft(pt,xt,𝜻t)) d𝛷t(𝜻t | xt) . (5) jS t(xt) ∶= N t ∑ k=1 𝛼k𝜑k(xt) , (6)  jS t(pt,x(k) t) ∶= u(ct(pt,x(k) t)) + 𝜌𝔼t [ jS t+1(ft(pt,x(k) t,𝜻t)) ]. 2 In this paper, we assume that numerical quadrature is possible. The applicability of the algorithms does not change if the distribution is generated by a different method, e.g., a Monte Carlo simulation. However, due to their slow convergence, Monte Carlo integration methods need a large simulated sample to obtain the accuracy required for the numerical optimization routine, see, e.g., the discussion by Cai (2019). 191 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… (Valentin 2019), the computational effort needed for the computation of the coefficients can be further reduced by using the unidirectional principle. This is facilitated by weakly fundamental not-a-knot splines and the insertion of some additional grid points. 3.1 Not‑A‑Knot B‑Spline Basis Let p∈ℕ0 , m∈ℕ , and 𝝃=(𝜉0,…,𝜉m+p) be an increasing sequence of real numbers. The B-spline bp k,𝝃 of degree p for the knot sequence 𝝃 is defined via the Cox–de Boor recurrence (Cox 1972; de Boor 1972; Höllig and Hörner 2013) where k=0, …,m−1 . It can be shown that for m>p , the B-splines bp 0,𝝃 ,…,b p m−1,𝝃 form a basis of the spline space Sp 𝝃 ∶= span{b p k,𝝃 ∣k=0, …,m−1 } on Dp 𝝃 ∶= [𝜉p,𝜉m ] (Höllig and Hörner 2013). The space Sp 𝝃 contains exactly those functions s ∶D p 𝝃 → ℝ which are piecewise polynomial of degree smaller or equal to p on every knot interval [𝜉k,𝜉k+1] in the interior of Dp 𝝃 ( k=p,…,m−1 ) and at least p−1 times continuously differentiable at every knot 𝜉k in the interior of Dp 𝝃 , k=p+1, …,m−1 (Höllig and Hörner 2013). For simplicity, we restrict the considerations and results in this paper to cubic B-splines (i.e., p=3 ) although it is important to note that our method can be generalized to arbitrary odd B-spline degrees. A common special case is the case of linear B-splines ( p=1 , so called hat functions), which are commonly used as basis functions for sparse grids. We consider equidistant grid points x 𝓁 ,i∶= ih𝓁 on the unit interval [0,1] where 𝓁∈ℕ0 is the level, i=0, …,2 𝓁 is the index, and h𝓁∶= 2−𝓁 is the mesh size. We want to find basis functions 𝜑 𝓁 ,i∶[0, 1] → ℝ such that we can interpolate a given objective function f∶[0, 1] → ℝ on the equidistant grid of level 𝓁 by a linear combination of the basis functions: with some 𝛼 𝓁 ,i∈ℝ . The most straightforward choice of B-splines for 𝜑 𝓁 ,i are uniform B-splines that are scaled and translated versions of the cardinal B-spline b3 ∶= b 3 0,(0,1,2,3,4) (Pflüger 2010; Valentin 2019): (7) b p k,𝝃(x) ∶= x−𝜉k 𝜉k+p−𝜉k bp−1 k,𝝃(x)+ 𝜉 k+p+1 −x 𝜉k+p+1−𝜉k+1 bp−1 k+1,𝝃(x) , b 0 k,𝝃(x) ∶= { 1 if x∈[𝜉k,𝜉k+1), 0 otherwise, (8)  f(x𝓁,j)=f(x𝓁,j)for all j=0, …,2 𝓁, where  f∶= 2 𝓁 ∑ i=0 𝛼𝓁,i𝜑𝓁, i (9) 𝜑unif 𝓁 ,i (x) ∶= b 3 (2 𝓁 x+2−i) . 192 P.Schober et al. 1 3 The resulting uniform B-splines 𝜑unif 𝓁 ,i , i=0, …,2 𝓁 , are exactly the B-splines b3 k,𝝃unif , k=0, …,m−1 , that arise from Eq. (7) when choosing the uniform knot sequence 𝝃unif ∶= (x𝓁 ,−2 ,x𝓁 ,−1 ,…,x𝓁 ,2 𝓁 +2) and m∶= 2𝓁+1 . However, the interpolation domain on which the B-splines span the spline space would only be D3 𝝃 unif =[𝜉p,𝜉m]=[x𝓁,1,x𝓁,2𝓁−1]=[2 −𝓁 ,1−2 −𝓁] . This interval does not contain the two boundary grid points x 𝓁 ,0 =0 and x 𝓁 ,2 𝓁 =1 . This leads to interpolation problems since the spline space on [0,1] is not contained in the spanned space of the basis functions 𝜑unif 𝓁 ,i . Even simple polynomials such as f(x)=4(x−0.5)2 cannot be represented exactly with the B-spline basis on the whole domain [0,1] as shown by Fig.1a and Valentin and Pflüger (2016). Consequently, the approximation quality for more complex functions like the value functions we interpolate in this paper deteriorates unnecessarily, which means that the economic results are not as conclusive as they could be. As a remedy, we impose so-called not-a-knot boundary conditions by removing the left-most and right-most inner grid points x 𝓁 ,1 and x 𝓁 ,2 𝓁 −1 from the knot sequence (Höllig and Hörner 2013; Valentin 2019). To keep the number m=2𝓁+1 of B-splines the same, we have to insert two additional knots outside the domain: The new interpolation domain D3 𝝃 nak =[x𝓁,0,x𝓁,2𝓁 ] is now the whole unit interval [0,1], containing all grid points at which we interpolate. As a result, the not-a-knot B-spline functions (10) 𝝃nak ∶= (x𝓁 ,−3 ,…,x𝓁 ,0 ,x𝓁 ,2 ,…,x𝓁 ,2 𝓁 −2 ,x𝓁 ,2 𝓁,…,x𝓁 ,2 𝓁 +3 ) . (11) 𝜑nak 𝓁,i∶= b 3 i,𝝃 nak ,i=0, …,2 𝓁, (a) (b) Fig. 1 Nodal B-spline bases and interpolation of parabola. a Nodal uniform B-spline basis of level 3 in 1D and interpolation of the parabola f(x)=4(x−0.5)2 with this basis, resulting in oscillations near the boundary. b The corresponding nodal not-a-knot B-spline basis interpolates the parabola f exactly. The knots at the left-most and right-most inner grid points x𝓁,1 and x𝓁,7 (crosses) are removed 199 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… interpolant (24) follows in Sect.4.4. Major parts from the Sects.4.1–4.4 are taken from our previous work in the recently submitted Ph.D. thesis of Valentin (2019). 4.1 Solution fortheValue Function Algorithm 1 shows solveValueFunction, generating the value function interpolants j S,1 t and j S,p t ( t=0, …,T ). The algorithm follows a simple optimize–refine–interpolate scheme, which is presented in Fig. 6: First, Eq. (23) is solved on an initial sparse grid (optimize). Then, we refine the grid spatially adaptively. Finally, the resulting grid data are interpolated with hierarchical higher-order B-splines. At the beginning of every iteration t, the grid of the piecewise linear interpolant is reset to an initial, possibly regular sparse grid. It would also be possible to reuse the grid from the previous iteration t+1 . However, the results we then obtain become worse, likely due to the different characteristics of j S,1 t for different t (e.g., kinks). 4.2 Optimization The optimize step can be seen in Algorithm2. This algorithm accepts in j S,1 t a spatially adaptive sparse grid 𝛺 S t ={x (k) t ∣k=1, …,N t} where the function values j S,1 t (x (k) t) may already be known for some grid points x(k) t if optimize is called from within refine. The function optimize computes the missing value function values. For t=T , we assume that the terminal solution jT can be computed by some Fig. 6 Scheme of the generation of value function interpolants with solveValueFunction (left, Algorithm 1), which repeatedly calls the optimize algorithm (right, Algorithm 2), which in turn consists of various sub-functions. The function optimize iterates over all state grid points x t =x (k) t ( k=1, …,Nt ) and calls optimizeSinglePoint for each point. The optimization method evaluates the objective function and its gradient at a sequence of different policy points pt to find p opt,S t (x ( k ) t) . This evaluation (denoted by evalObjFcnGrad) has to implicitly compute the expectation in Eq.(23), which is done using a quadrature rule. For every quadrature point 𝜻t =𝜻 (j) t ( j=1, …,Qt ), evalQuadPoint computes the corresponding value of the expression in the expectation. Finally, evalInterpPartDeriv evaluates the interpolant j S,p t+1 and its partial derivatives for which we have to loop over the state dimensions o=1, …,d 200 P.Schober et al. 1 3 function computeKnownTerminalSolution.3 Otherwise, for t<T , we solve the maximization problem(23) by using the higher-order B-spline interpolant j S,p t+1 of the previous iteration t+1 (optimizeSinglePoint). The computations for the different x(k) t are independent of each other, which means that they can be computed in parallel (Cai etal. 2015; Horneff etal. 2016).4 After generating all missing data, we update the hierarchical surpluses of the piecewise linear interpolant j S,1 t to interpolate the new data at all grid points of 𝛺S t . 4.3 Refinement For adaptive refinement, the criterion is the common surplus-volume (Pflüger 2010). We use the piecewise linear interpolant for the surplus-based grid generation as the surpluses are easier to compute in the piecewise linear case, and as they are more meaningful due to the integral representation formula (Bungartz and Griebel 4 Such a problem is usually referred to as embarrassingly parallel. 3 In any case, the terminal solution may be computed as the solution of the corresponding single-time optimization problem, e.g., jT(x (k) T )=maxp T {u(cT(x (k) T ,pT)) } . 201 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… 2004). Algorithm 3 shows how to generate the spatially adaptive sparse grid in solveValueFunction (Algorithm1). Parameters are the tolerance 𝜀∈ℝ≥0 by which the set of grid points to be refined is determined and the number q∈ℕ0 of refinement iterations. 4.4 Solution fortheOptimal Policies To construct the optimal policies, we use the higher-order B-spline interpolant jS,p t and the optimal policies popt t(x(k) t) at the grid points x(k) t ( k=1, …,Nt ) obtained from Algorithm1. We then spatially adaptively refine the grid for each policy to construct a policy interpolant of degree one, popt,S,1 t , for each t=1, …,T . The corresponding Algorithm1 is similar to solveValueFunction (Algorithm4), except that it operates on the policy interpolants instead of the value function interpolant. The functions optimize and refine have been replaced by corresponding policy versions optimizePolicy and refinePolicy that work very much like their value function counterpart. In the optimization step, optimizePolicy only has to generate new values if the initial regular sparse grid for the policies is not contained in the grid of jS,p t . The policy grid is then refined independently of the value function grid. The iterations are independent of each other, which means that they can be parallelized.5 5 In principle, one could generate policy interpolants of degree p, popt,S,p t , by adding an extra interpolation step after refinePolicy (Valentin 2019). 202 P.Schober et al. 1 3 5 Application: Transaction Costs Problem First, we introduce the dynamic portfolio choice model with transaction costs (Sect.5.1). We then describe its numerical solution in Sect.5.2 and derive our error measure (unit-free Euler equation errors) in Sect.5.3. We verify our solution economically on a two-dimensional full grid (Sect.5.4). Then, we analyze the time complexity of the solution approach and show the impact of the choice of basis functions on the computational complexity in Sect.5.5. Therefore, we solve the problem with B-splines of cubic degree on a regular sparse grid and compare them to the linear approach as used by Brumm and Scheidegger (2017). We find that we save approximately one order of magnitude in computational complexity with cubic B-splines already for three dimensions. Furthermore, we compare the results of our approach based on analytical gradients with the results we obtain with finite differences. For d>3 , solutions on regular sparse grids are no longer feasible to compute in suitable numerical accuracy. In Sect.5.6, we illustrate how spatial adaptivity allows us to solve the transaction costs problem up to d=5 accurately by showing pointwise error decay and convergence of our approach. Finally, we present in Sect.5.7 economic results for the transaction costs problem in higher dimensions. The results in these sections are also contained in the recently submitted Ph.D. thesis of Valentin (2019). 5.1 Transaction Costs Problem As the transaction costs problem is easiest described in vector notation, let us denote the unit vector with 1 . For two vectors a , b , we define the Hadamard product as (a⊙b)i∶= aibi for all i. In the transaction costs problem the investor maximizes expected utility from consumption(1). Therefore, at time t, she tracks her wealth Wt∈ℝ≥0 and fractions of wealth xt∈[0, 1]d invested in stocks. Her choices are how much to buy of the d stocks 𝜟+ t ∈ℝ d ≥0 with transaction costs 𝜏 𝜟 + t or sell 𝜟− t ∈ℝ d ≥0 with transaction costs 𝜏𝜟− t where 𝜏>0 is a cost factor. Additionally, she can invest in a transaction-costfree money market account Bt , yielding a risk-free return rf∈ℝ≥0 . We assume the returns rt∈ℝ≥0 that the d stocks earn from t to t+1 are independent and identically 203 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… lognormally distributed with mean 𝝁 and covariance matrix 𝜮 : rt∼LN(𝝁,𝜮) (Cai 2009; Cai and Judd 2010). The investor’s consumption Ct in period t is the residual of her wealth that is not invested in stocks or bonds, reduced by the transaction costs for rearranging her portfolio in this period: The state dynamics from t to t+1 are thus given by: The investor faces the optimization problem with utility function u from Eq.(2) subject to the constraint for all t=0, …,T , where 𝜟+ t≥0 , 𝜟− t∈[0,xtWt] , Bt≥0 , and 1⊤ ⋅ xt≤1 . Here, a minimum consumption level Cmin must be maintained, and the final stock holdings xTWT are assumed to be sold before they can be consumed. In addition, at no point in time t the investor can sell more of the stocks than her current holdings xtWt . The problem can be simplified by normalizing the value function jt=Jt∕Wt , consumption ct=Ct∕Wt , and investment choices bt=Bt∕Wt , 𝜹+ t =𝜟 + t ∕W t , 𝜹− t=𝜟− t∕Wt with respect to wealth Wt for each t. The investor’s normalized consumption ct in period t is then The state dynamics from t to t+1 can be expressed in terms of the portfolio value in t+1 : (25) Ct = ( 1−1 ⊤ ⋅x t) W t −B t −(1+𝜏)1 ⊤ ⋅𝜟 + t −(𝜏−1)1 ⊤ ⋅𝜟 − t. (26a) Wt+1 =B t r f + ( x t W t +𝜟+ t −𝜟− t)⊤ ⋅r t, (26b) x t+1= ( xtWt+𝜟 + t−𝜟 − t ) ⊙rt Wt+1 . (27a) Jt(Wt,xt)= max B t ,𝜟+ t ,𝜟− t { u(Ct)+𝜌𝔼t [ Jt+1 ( Wt+1,xt+1 )]} ,t<T , (27b) J T (W T ,x T )=u (( 1− 𝜏 1 ⊤ ⋅x T) W T), (27c) Bt +(1+𝜏)1 ⊤ ⋅𝜟 + t +(𝜏−1)1 ⊤ ⋅𝜟 − t≤( 1−1 ⊤ ⋅x t) W t −C min , (28) ct =1−1 ⊤ ⋅x t −b t −(1+ 𝜏 )1 ⊤ ⋅𝜹 + t −( 𝜏 −1)1 ⊤ ⋅𝜹 − t. (29) 𝜋t+1 ∶= b t r f +(x t +𝜹 + t −𝜹 − t ) ⊤ ⋅r t (30a) Wt+1=Wt𝜋t+1, (30b) xt+1= ( xt+𝜹 + t−𝜹 − t ) ⊙rt 𝜋t+1 . 204 P.Schober et al. 1 3 With this normalization, the solution to problem (27a–27c) can be expressed as Jt =W 1−𝛾 t j t for each t=0, …,T with subject to the constraints for all t=0, …,T , where the minimum consumption level cmin =Cmin∕Wt is also normalized with respect to wealth (see "Appendix A.2"). Now the investor’s optimization problem no longer depends on Wt , and hence one state variable can be eliminated. The nonnormalized optimal choices can be obtained by multiplication with a given wealth Wt for any t and state xt . 5.2 Numerical Solution To compute the solution to the transaction costs problem, we use the certainty equivalent transformation  j t of the normalized value function jt , which reduces the curvature of the value function when the utility is of Constant Relative Risk Aversion type, Eq. (2) (Garlappi and Skoulakis 2009). Since this transform is strictly monotone, any maximizer of  jt also maximizes jt . The optimization problem then reads (31a) jt(xt)= max b t, 𝜹+ t, 𝜹− t{ u(ct)+𝜌𝔼t [ 𝜋1−𝛾 t+1jt+1(xt+1) ]} ,t<T , (31b) j T (x T )=u ( 1− 𝜏 1 ⊤ ⋅x T), (31c) bt +(1+𝜏)1 ⊤ ⋅𝜹 + t +(𝜏−1)1 ⊤ ⋅𝜹 − t≤ 1−1 ⊤ ⋅x t −c min , (31d) 𝜹+ t≥ 0 , (31e) 𝜹− t≥0, (31f) 𝜹− t≤xt, (31g) bt≥0, (31h) 1⊤ ⋅ xt≤1, (32)  j t (x t )= ( (1− 𝛾 )j t (x t ) ) 1 1−𝛾 , (33a)  jt(xt)= max b t ,𝜹+ t ,𝜹− t{( c1−𝛾 t+𝜌𝔼t [( 𝜋t+1 jt+1(xt+1) ) 1−𝛾 ]) 1 1−𝛾 }, 205 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… with constraints Eq.(31c) to (31h) (see "AppendixA.3"). The optimization problem for a given state xt is solved with the SQP solver SNOPT (Gill et al. 2005), see "Appendix A.4" for the specific objective function and its gradient. Since the distribution of the returns rt is multivariate lognormal and state-independent, we can compute the expectation in Eq.(6) using Gauss–Hermite quadrature. For this, we also use a sparse grid quadrature rule, thus breaking the curse of dimensionality when including stochastic risk factors (see the appendix of Horneff etal. 2016 for details). The constraint (31h) constrains the state space. The resulting eligible subspace 𝛺Simplex ={x t ∈[0, 1]d∣1 ⊤ ⋅x t≤ 1}⊂[0, 1]d is a d-dimensional simplex, not a rectangular domain as needed for the sparse grid approximation. We solve this problem by assuming that any state attained that is not eligible is cropped to an eligible state by selling all stock holdings pro rata until all constraints are satisfied. That is, money is transferred from stocks to wealth, for which the proportionate transaction costs are deducted (see "AppendixA.5"). The approximation of the value function is then evaluated at this eligible state. The optimization ran over an investment horizon of T=6 years and had a fixed period length of 1year. The risk aversion 𝛾=3.5 , the risk-free rate rf=4% , and the transaction costs factor 𝜏=1% were taken from Cai and Judd (2010). We extended the return distribution parametrization of Cai and Judd (2010) to five dimensions: and set the time discount factor to 𝜌=0.97 and cmin =0.001 to ensure that minimal consumption was taking place. For d stocks we used the first d entries of 𝝁 and the elements 𝜮i,j , i,j≤d , as the return distribution parametrization. As initial grids, we used regular sparse grids 𝛺S n,d of level n and dimension d. The code was written in MATLAB where the interpolation on sparse and full grids was implemented by a MEX file interface to the sparse grids C++ toolbox SG++ (sgpp.spars egrid s.org, Pflüger 2010). The quadrature routine was implemented by a MEX file interface to the TASMANIAN sparse grids C++ toolbox (Stoyanov 2017) as TASMANIAN allows us to integrate a real valued function over a Gaussian density using Hermite polynomials on sparse grids (see the appendix of Horneff etal. 2016). We used the SNOPT implementation of the Numerical Algorithms Group (www.nag.co.uk). If convergence of the optimizer was not observed, we stopped the optimization after 100 iterations. To avoid being stuck in local minima, we repeated the optimization process for a varying number of initial multi-start points (in the range of a few dozens). All computations were performed on the compute cluster LOEWE-CSC (csc.uni-frank furt.de) where we exclusively allocated three compute nodes with two Intel Xeon (33b)  jT(xT)=1−𝜏 1 ⊤ ⋅ xT, (34) 206 P.Schober et al. 1 3 E5-2670 v2 CPUs (ten cores at 2.5 GHz, 20 threads) each, i.e., 120 threads in total, and 4000 MB RAM per thread. 5.3 Error Measurement Any optimal policy p opt t ∶= (b opt t ,𝜹 +,opt t ,𝜹 −,opt t ) ⊤ must satisfy the first order conditions of the Lagrangian at any given state xt for each t<T . Specifically, for the transaction costs problem—when neglecting binding constraints—we obtain from the first-order condition with regard to the optimal bond policy bopt t : where Rearranging Eq.(35) and taking the root ( ⋅ )−1∕𝛾 yields the unit-free Euler equation error, which should be 0 for any given state xt in the eligible domain 𝛺Simplex . However, the state space cropping distorts unit-free Euler equation errors. This is due to three sources: Firstly, the cropping already occurs for large stock holdings 1⊤ ⋅ xt that are less than one as stocks have to be sold to maintain minimum consumption cmin . Secondly, transaction costs for selling the stocks are deducted. Thirdly, even if neither minimum consumption is required, nor transaction costs are incurred the error at the hyperplane 1⊤ ⋅ xt=1 does not vanish even for full grid solutions. Only in the limit, as the resolution of the grid goes to infinity, the error will vanish. Economically, the region near this hyperplane is not significant as such large stock fractions are unusual, which is confirmed by Monte Carlo simulations. We therefore use the weighted Euler equation error (35) − copt t −𝛾+𝜌𝔼t [( 𝜋opt t+1 jS t+1 ) −𝛾rf (  jS t+1− ( 𝛁xt+1  jS t+1 ) ⊤ ⋅xopt t+1 )] = 0, (36a) copt t =1−1 ⊤ ⋅x t −b opt t −(1+ 𝜏 )1 ⊤ ⋅𝜹 +,opt t −( 𝜏 −1)1 ⊤ ⋅𝜹 −,opt t, (36b) 𝜋opt t+1 =b opt t r f +(x t +𝜹 +,opt t −𝜹 −,opt t ) ⊤ ⋅r t, (36c) xopt t+1= ( xt+𝜹 +,opt t−𝜹 −,opt t ) ⊙rt 𝜋opt t+1 . (37) 𝜀 t(xt)= ( 𝜌𝔼t [( 𝜋opt t+1 jS t+1 ) −𝛾rf (  jS t+1− ( 𝛁xt+1  jS t+1 ) ⊤ ⋅xopt t+1 ) copt t 𝛾 ]) − 1 𝛾 − 1, 207 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… instead of 𝜀t . Alternatives would be restricting the state domain in which the error is computed or weighting the error with the probability that a given state occurs in a Monte Carlo simulation.6 We then choose the same N=1000 points x(k) ∈𝛺 Simplex ( k=1, …,N ) for all times t=0, …,T−1 and compute the errors 𝜀w t (x (k)) for each t.7 We report the L2 norm scaled by √d! and the L∞ norm for each t: For details on the error derivation see "AppendixA.6". In principle, we could also compare the solutions  j S t and optimal policies p opt,S t obtained on sparse grids with the full grid solution, e.g., in a point-wise way. However, full grid solutions with acceptable resolutions are computationally infeasible already in d>2 . In addition, the Euler equation error does not compare numerical solutions with each other, but rather measures the accuracy of any solution, regardless of whether it is obtained numerically or analytically. 5.4 Economical Verification We show in Fig. 7 a full grid solution for the case of d=2 stocks, i.e., { x (k) t ∣k=1, …,N t }={0, 2−n,…,1}d for some fixed level n∈ℕ (here, n=7 and Nt =(2 7 +1) 2 = 16641 ) and for all t=0, …,T . The red dot (xt,1,xt,2)=(0.1509, 0.1831) shows the so-called Merton point for which Merton (1969) derives that in the case of 𝜏=0 the optimal stock fractions xopt t are constant over time and wealth. When faced with transaction costs, Magill and Constantinides (1976) find that the investor must weigh up the benefits of improved diversification against the associated transaction costs for rebalancing the portfolio. This leads to the no-trade region (red outline). If xopt t lies within this region, the investor does not alter her portfolio. In discrete time consumption and portfolio choice, the no-trade region is known to be a convex set, and, if the current stock fraction is outside this region, the optimal policy is to move to the convex hull (38) 𝜀w t (x t ) ∶= ( 1−1 ⊤ ⋅x t) 𝜀 t (x t) (39a) 𝜀 w,L2 t∶= √ √ √ √ 1 N N ∑ k=1 | 𝜀w t(x(k)) | 2 , (39b) 𝜀w,L∞ t ∶= max{ | 𝜀 w t (x (k) ) | ∣k=1, …,N} . (40) xopt t∶= 𝜮 −1 (𝝁−1rf) 𝛾 , 6 Another possibility is to transform the non-rectangular domain to the unit hypercube, e.g., by analyzing the principal components of the ergodic distribution (Judd etal. 2014). 7 We choose 1000 ×d! points from [0, 1]d and discard all points that are not in 𝛺Simplex . Here, the d! corrects for the volume of the d-dimensional simplex obtained by cropping. Thus, there are only N≈1000 points in 𝛺Simplex . 208 P.Schober et al. 1 3 of the set (Abrams and Karmarkar 1980; Constantinides 1979). Naturally, the Merton point lies inside the no-trade region. We can confirm this result for our optimal policies, i.e., if we choose any point outside the no-trade region (but within the eligible domain 𝛺Simplex ), computing xt+ 𝜹 +,opt t− 𝜹 −,opt t then results in a boundary point of the no-trade region. Figure7 also shows the impact of the state space cropping as the eligible subspace 𝛺Simplex ⊂[0, 1] d is not a rectangular domain as needed for the sparse grid interpolation. This is the reason why the certainty equivalent value function  j S t is zero in [ 0, 1] d ⧵𝛺 Simplex and the optimal sell policies 𝜹−,opt,S t contain a diagonal kink at the hyperplane 1⊤ ⋅x t . Obviously, computing full grid solutions is only computationally feasible for low dimensionalities d due to the curse of dimensionality. The two-dimensional solution of level n=7 took over nine hours to compute on the LOEWE-CSC cluster with 120 threads. The solution of the next level is estimated to already take one week. Hence, full grid solutions can only be computed up to d=3 due to prohibitively long computational times for d≥4 . This underlines the need for sophisticated discretization techniques such as sparse grids. 5.5 Savings inComplexity Using B‑Splines A complexity analysis reveals that the difficulty of solving dynamic portfolio choice models quickly grows with the dimensionality d: The number of necessary arithmetic operations grows like (see Fig.6) Fig. 7 Full grid solution for the transaction costs problem with d=2 stocks. Shown are the certainty equivalent value function  j S t (x t) (top left) and the optimal policies popt,S t(xt) for t=0 . Also shown are the Merton point xopt t (red dot) and the no-trade region (red outline) 215 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… Finally, we present in Table2 the computational times and numerical errors of the sparse grid solutions underlying Fig.13 and Table1. (a) (b) (c) (d) Fig. 13 Average values of wealth Wt (blue), unnormalized optimal bonds Bt (purple), unnormalized optimal consumption Ct (green), and unnormalized stock holdings  xt ∶= (x t +𝜹 + t −𝜹 − t )W t after buying and selling (red) in a Monte Carlo simulation of 105 individuals where we assume that W0=$1 for all individuals. In addition, the plots show the evolution of the scaled L2 error 𝜀 w,L 2 t over time t (gray, right axes) Table 1 Simulated stock fractions x t,o ∕(1 ⊤ ⋅ x t) for t=0 and stock fractions implied by the Merton points xopt t,o ∕(1 ⊤ ⋅x opt t) ( o=1, …,d ) for the Monte Carlo simulations obtained by evaluating the optimal policy interpolants computed on spatially adaptive sparse grids d xt,o ∕(1 ⊤ ⋅ x t) xopt t , o∕(1⊤ ⋅x opt t ) 2 (0.441,0.559) (0.452,0.548) 3 (0.300,0.317,0.383) (0.314,0.302,0.384) 4 (0.239,0.238,0.253,0.270) (0.275,0.185,0.250,0.289) 5 (0.199,0.188,0.197,0.205,0.212) (0.275,0.122,0.176,0.203,0.223) 216 P.Schober et al. 1 3 6 Conclusion In this paper, we are the first to develop an approach to accurately solve high-dimensional dynamic portfolio choice models in discrete time that require smooth approximations or gradient-based optimization. With our approach, we have addressed the three key issues of solving these models by means of value function iteration: the curse of dimensionality, the lack of spatial adaptivity, and the lack of continuous gradients all at once by using B-splines on sparse grids with spatially adaptive refinement. We have solved a dynamic portfolio and consumption choice model with transaction costs to study the numerical accuracy of our approach. Solutions to the transaction costs problem with value function iteration have achieved economically acceptable results already for lower resolutions of the interpolation grid than in our presented example. Our approach, however, can easily be applied to other dynamic portfolio choice models or any high-dimensional economic model that require such a high resolution. We have solved the transaction costs problem with up to five stocks and one risk-free bond, i.e., a five-dimensional interpolation and an eleven-dimensional optimization problem per time step. Using spatially adaptive refinement of the optimal policies, we have obtained maximum unit-free Euler equation errors around 5% for the five-dimensional problem and even lower maximum errors for lower-dimensional problems. This showcases the high accuracy of the proposed spatially adaptive solution scheme for the optimization of continuous choices, which relies on smooth approximations of the value function and the gradient. We have shown convergence of our approach in up to four dimensions. Here, spatially adaptive refinement of the optimal policies decreased the maximum Table 2 Number of stocks d, base level n, refine tolerance 𝜀 , grid points of the base grid | 𝛺 S n,d| for d stocks and level n, grid points Nt of the refined grid, added points 𝛥Nt , computational time and weighted Euler equation errors 𝜀 w,L 2 t , 𝜀w,L∞ t , at t=0 . For the optimal policy p opt,S,1 t rows, the number Nt is the average number 1 ∕mp ∑m p j=1 Nt, j of grid points over all policy grids for t=0 where Nt,j is the number of grid points of the j-th policy entry dInterp. n 𝜀 | 𝛺 S n,d| Nt 𝛥 N t Time 𝜀 w,L 2 t 𝜀w,L∞ t 2  j S,1 t 4 ∞ 113 113 0 0 min 8.9 e-06 9.0 e-05 p opt,S,1 t 43 e-09 113 5154 5041 3 min 3 j S,1 t 4 ∞ 593 593 0 3 min 3.7 e-04 3.6 e-03 p opt,S,1 t 42 e-07 593 28735 28142 5 h 19 min 4  j S,1 t 4 ∞ 2769 2769 0 3 h 41 min 1.5 e-02 3.9 e-02 p opt,S,1 t 42 e-04 2769 3343 574 3 h 30 min 5 jS,1 t 4 ∞ 12033 12033 0 9 h 23 min 2.0 e-02 4.8 e-02 popt,S,1 t 43 e-04 12033 12572 539 3 h 5 min 217 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… Euler equation error by nearly two orders of magnitude in the four-dimensional case compared to regular sparse grids without spatially adaptive refinement. This has shown that only with spatial adaptivity, high-dimensional problems can be solved accurately. Finally, we have given a rigorous analysis of the complexity of our approach for dynamic portfolio choice models in general, not only for the transaction costs problem for which we have verified our analysis of complexity with measurements of computational time. We have found that the sole availability of the gradient for the optimization process has saved nearly one order of magnitude in computational complexity in the three-stock case. We expect even larger reductions of the computational complexity for higher-dimensional problems. Compared to finite differences with interpolation on hat functions as used by Brumm and Scheidegger (2017), we have saved considerably more than one order of magnitude in computational complexity and one order of magnitude in total computational time in three dimensions. There are certain limitations to the applicability of spatially adaptive sparse grids to solve high-dimensional dynamic economic models: Firstly, sparse grid approximations are not shape-preserving, which is especially of importance for value function iteration with interpolation (Cai and Judd 2012). Secondly, the calculation of the coefficients of the B-spline interpolant is time-consuming and not trivial to parallelize since the solution to a system of linear equations has to be computed in every time step. Thirdly, the exact choice of the refinement tolerance for value function and policy interpolants is subject to trial-and-error. Choosing a refinement tolerance that is too low will lead to too many points that are inserted and may cause instability of the entire scheme if the optimizer does not give perfect results. Future improvements of our approach may lie in the use of problem-tailored adaptivity criteria (Brumm and Scheidegger 2017; Pflüger 2012) instead of the simple surplus-based refinement criterion. Appendix A A.1 Sequential Quadratic Programming (SQP) SQP methods are well-suited to compute the solution to problem(3). These methods use the linearization of the Lagrangian with Kuhn-Tucker multipliers 𝝀t∈ ℝ mg to set up a quadratic programming problem (44) Lt (p t ,𝝀 t ,x (k) t ) ∶=  jS t (p t ,x (k) t )+𝝀⊤ t ⋅g t (p t ,x (k) t ) , (45) max d (i) t { 𝛁pt  jS t(p(i) t,x(k) t)⊤ ⋅d(i) t+ 1 2d(i) t ⊤ ⋅𝛁2 pt  jS t(p(i) t,x(k) t)⋅d(i) t }, 218 P.Schober et al. 1 3 which finds the search direction d(i) t for the current iterate p (i) t starting from an initial guess p (0) t ∈ 𝛹⊂ ℝm p .8 This problem can hence be solved by means of standard quadratic programming where frequently the Hessian 𝛁2 pt  j S t is approximated by the BFGS method (Fletcher 2013). Finally, the next iterate p (i+1) t is chosen by an appropriate line search procedure over the step length 𝜎i : If the gradient 𝛁pt  j S t is not available, implementations of SQP methods approximate it by finite differences. This leads to mp additional evaluations of  j S t per quadratic programming iteration and to 2mp additional evaluations if central finite differences are used. Some SQP routines allow the user to choose finite difference approximations, some automatically use central finite differences in certain situations to achieve higher accuracy when needed. However, almost all SQP routines allow the user to provide the gradient 𝛁pt  j S t (and many also the gradients of the constraints) to save computing time and to increase accuracy. For details see Fletcher (2013), Gill etal. (2005). The gradient 𝛁pt  j S t of the target function(6) at the grid point x (k) t , can be supplied to the SQP routine by evaluating the approximation of the gradient 𝛁 f t j S t+1 rather than using finite differences. A.2 Proof oftheNormalization Theorem1 For all t=0, …,T, where Jt is the solution of problem (27a–27c) and jt the solution of problem (31a–31h). Proof We have u( C t) =u ( W t c t) =W 1−𝛾 t u ( c t) for all t due to the choice of the utility function(2). Base case t=T : Because of Eqs.(27b) and (31b) it is Inductive hypothesis: For t+1 it is Jt+1( W t+1 ,x t+1) =W 1−𝛾 t+1 j t+1 (x t+1) . (46) p (i+1) t ∶= p (i) t + 𝜎i d (i) t. (47) 𝛁 pt  jS t=𝛁ptu(ct(pt,x(k) t)) + 𝜌𝔼t [ 𝛁ftjS t+1(ft(pt,x(k) t,𝜻t))⊤⋅𝛁ptft(pt,x(k) t,𝜻t) ], (48) J t (W t ,x t )=W 1−𝛾 t j t (x t ) , (49) J T(WT,xT)=u (( 1−𝜏1⊤⋅xT ) WT ) =W 1− 𝛾 Tu ( 1−𝜏1⊤⋅xT ) =W1−𝛾 T j T (x t ). 8 We define 𝛁xf∶= (𝜕fj∕𝜕xi)i,j (i.e., the transposed Jacobian) and “ ⋅ ” denotes the matrix-vector product. 219 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… Inductive step t+1→t : By the inductive hypothesis and Eq.(30a) it is ◻ A.3 Proof oftheCertainty Equivalent Transformation Theorem2 For all t=0, …,T, where jt the solution of problem (31a–31h) and  jt the solution of problem (33a, 33b). Proof Base case t=T : Because of Eqs.(31b) and (33b) it is Inductive hypothesis: For t+1 it is j t+1 (x t+1 )=1∕(1− 𝛾 )  j t+1 (x t+1 )1 − 𝛾. Inductive step t+1 → t : By the inductive hypothesis and Eq.(33a) it is where we used that ( ⋅ ) 1 ∕( 1 − 𝛾 ) is a strictly monotonously decreasing function and 1−𝛾<0 . ◻ (50) J t(Wt,xt)= max Bt,𝜟+ t,𝜟− t { u(Ct)+𝜌𝔼t [ W1−𝛾 t+1jt+1(xt+1) ]} =W1−𝛾 tmax bt,𝜹+ t,𝜹− t{u(ct)+𝜌𝔼t[𝜋1−𝛾 t+1jt+1(xt+1)] } =W 1−𝛾 tjt(xt). (51) jt(xt)= 1 1−𝛾  jt(xt)1−𝛾 , (52) jT(xT)=u ( 1−𝜏1⊤⋅xT ) = 1 1−𝛾 ( 1−𝜏1⊤⋅xT ) 1− 𝛾 =1 1−𝛾  jT(xt)1−𝛾. (53) jt(xt)= max bt,𝜹+ t,𝜹− t { u(ct)+𝜌𝔼t [ 𝜋1−𝛾 t+1jt+1(xt+1) ]} =max bt,𝜹+ t,𝜹− t{1 1−𝛾c1−𝛾 t+𝜌𝔼t[𝜋1−𝛾 t+1 1 1−𝛾  jt+1(xt+1)1−𝛾]} =1 1−𝛾min bt,𝜹+ t,𝜹− t{c1−𝛾 t+𝜌𝔼t[(𝜋t+1 jt+1(xt+1))1−𝛾]} =1 1−𝛾(max bt,𝜹+ t,𝜹− t{(c1−𝛾 t+𝜌𝔼t[(𝜋t+1 jt+1(xt+1))1−𝛾]) 1 1−𝛾})1−𝛾 =1 1−𝛾  jt(xt)1−𝛾, 220 P.Schober et al. 1 3 A.4 Analytical Gradients For t<T at state xt , let us define the objective function (6) for the certainty equivalent formulation of the transaction costs problem (33a, 33b) by where  j S t denotes the sparse grid B-spline approximation of Eq.(33a) evaluated at xt . With the sparse grid B-spline approximation of the gradient with respect to xt evaluated at xt , 𝛁xt  j S t , the gradient of the objective function (47) with respect to the policies bt , 𝜹+ t , and 𝜹− t is: A.5 State Space Cropping To obtain function values outside the feasible state space we virtually sell, if 1⊤ ⋅ xt>1 as many stocks as needed to meet the constraint 1⊤ ⋅ xt≤1 . We already might need to sell stocks even if 1⊤ ⋅ xt is smaller but close to one in order to satisfy the minimum consumption requirement (31c). In detail, we replace xt by  𝛽xt whenever  𝛽<1 where  𝛽>0 is a cropping factor that is determined by Here, ( 1 ⊤ ⋅x t −1 ⊤ ⋅(  𝛽x t ) ) is the amount of virtually sold stocks. Hence, the term in square brackets is the fraction of wealth that is still available after deducting the induced transaction costs. The product of this term with ( 1−1 ⊤ ⋅(  𝛽x t ) ) is the fraction of wealth that can be consumed after the virtual selling, which needs to be at least cmin . Solving Eq.(56) for  𝛽 and choosing the positive solution, we finally obtain (54)   jS t(bt,𝜹+ t,𝜹− t,xt) ∶= { c1−𝛾 t+𝜌𝔼t [( 𝜋t+1 jS t+1 ) 1−𝛾 ]} 1 1−𝛾 , (55a) 𝛁 bt   jS t=  jS t 𝛾 ( −c−𝛾 t+𝜌𝔼t [( 𝜋t+1 jS t+1 ) −𝛾rf (  jS t+1− ( 𝛁xt+1  jS t+1 ) ⊤ ⋅xt+1 )]), (55b) 𝛁 𝜹+ t   jS t=  jS t 𝛾 ( −(1+𝜏)1c−𝛾 t+𝜌𝔼t [( 𝜋t+1 jS t+1 ) −𝛾 rt⊙ ( 1 jS t+1+𝛁xt+1  jS t+1−1 ( 𝛁xt+1  jS t+1 ) ⊤ ⋅xt+1 )]), (55c) 𝛁 𝜹− t   jS t=  jS t 𝛾 ( −(𝜏−1)1c−𝛾 t−𝜌𝔼t [( 𝜋t+1 jS t+1 ) −𝛾 rt⊙ ( 1 jS t+1+𝛁xt+1  jS t+1−1 ( 𝛁xt+1  jS t+1 ) ⊤ ⋅xt+1 )]). (56) [ 1−𝜏 ( 1⊤⋅xt−1⊤⋅( 𝛽xt) )] ⋅ ( 1−1⊤⋅( 𝛽xt) ) =cmin . (57)  𝛽 = 𝜏 ( 1+1⊤⋅xt ) −1+ √ 𝜏2 ( 1−1⊤⋅xt ) 2 −2𝜏 ( 2cmin −1+1⊤⋅xt ) +1 2𝜏 1⊤⋅x t . 221 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… A.6 Euler Equation Error Derivation The Lagrangian (44) for the transaction costs problem in certainty equivalent formulation (33a, 33b) is given by: with   j S t from Eq.(54), 𝝀∈ℝ3d+2 and g t (b t ,𝜹 + t ,𝜹 − t ,x t )=(c t −c min ,𝜹 + t ,𝜹 − t ,x t −𝜹 − t ,b t ) ⊤ from constraints Eqs.(31c) to (31g). The first order condition with regard to bt is We neglect binding constraints, i.e., we assume 𝜆1 t =𝜆 3d+2 t = 0 , and set the error 𝜀w t(xt)=NaN whenever 𝜆1 t≠0 or 𝜆3d+2 t≠0 for any xt . Thus, we take the NaN -mean in Eq.(39a) and the NaN -maximum in Eq.(39b). Assuming 𝜆1 t =𝜆 3d+2 t = 0 in Eq.(59), plugging in Eq.(55a) for 𝛁bt   j S t at the optimum ( b opt t ,𝜹 +,opt t ,𝜹 −opt t ) ⊤ and dividing both sides by   jS t 𝛾 yields Eq.(35). Acknowledgements We thank the initiative High Performance Computing in Hessen for granting us computing time at the LOEWE-CSC cluster and the Lichtenberg High Performance Computer. We appreciate valuable remarks to improve this paper from two anonymous referees. Helpful comments have been provided by Johannes Brumm, Yannick Dillschneider, Kenneth Judd, and Alexander Ludwig. We thank Stefan Zimmer for helping us develop the weakly fundamental not-a-knot splines. Finally, Peter Schober thanks Raimond Maurer for supporting this research in every way possible. Funding Open Access funding enabled and organized by Projekt DEAL. This work was supported by funding from the German Investment and Asset Management Association (BVI), the Juniorprofessurenprogramm of the Landesstiftung Baden-Württemberg, and the DFG (Cluster of Excellence SimTech EXC310/EXC2075). Availability of data and materials/Code availability The code used to generate the results presented in this paper and all input data is publicly available at trix0r/BBSG. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/. (58) Lt (b t ,𝜹+ t ,𝜹− t ,𝝀 t ,x t ) ∶=   jS t (b t ,𝜹+ t ,𝜹− t ,x t )+𝝀⊤ t ⋅g t (b t ,𝜹+ t ,𝜹− t ,x t ) , (59) 𝛁bt L t =𝛁 bt   jS t −𝜆1 t 𝛁 bt copt t −𝜆3d+2 t = 0. 222 P.Schober et al. 1 3 References Abrams, R. A., & Karmarkar, U. S. (1980). Optimal multiperiod investment-consumption policies. Econometrica, 48(2), 333–353. Barberis, N., & Huang, M. (2009). Preferences with frames: A new utility specification that allows for the framing of risks. Journal of Economic Dynamics and Control, 33(8), 1555–1576. Barthelmann, V., Novak, E., & Ritter, K. (2000). High dimensional polynomial interpolation on sparse grids. Advances in Computational Mathematics, 12(4), 273–288. Bellman, R. (1954). The theory of dynamic programming. Bulletin of the American Mathematical Society, 60(6), 503–515. Bellman, R. (1961). Adaptive control processes: A guided tour. Princeton: Princeton University Press. https ://doi.org/10.1515/97814 00874 668. Brumm, J., & Grill, M. (2014). Computing equilibria in dynamic models with occasionally binding constraints. Journal of Economic Dynamics and Control, 38, 142–160. Brumm, J., & Scheidegger, S. (2017). Using adaptive sparse grids to solve high-dimensional dynamic models. Econometrica, 85(5), 1575–1612. Bungartz, H.-J., & Griebel, M. (2004). Sparse grids. Acta Numerica, 13, 147–269. Cai, Y. (2009). Dynamic programming and its application in economics and finance. Ph.D. thesis, Stanford University. https ://purl.stanf ord.edu/zd335 yg688 4. Cai, Y. (2019). Computational methods in environmental and resource economics. Annual Review of Resource Economics, 11(1), 59–82. Cai, Y., & Judd, K. L. (2010). Stable and efficient computational methods for dynamic programming. Journal of the European Economic Association, 8(2–3), 626–634. Cai, Y., & Judd, K. L. (2012). Shape-preserving dynamic programming. Mathematical Methods of Operations Research, 77(3), 407–421. Cai, Y., & Judd, K. L. (2015). Dynamic programming with Hermite approximation. Mathematical Methods of Operations Research, 81(3), 245–267. Cai, Y., Judd, K. L., Thain, G., & Wright, S. J. (2015). Solving dynamic programming problems on a computational grid. Computational Economics, 45(2), 261–284. Cai, Y., Judd, K. L., & Xu, R. (2020). Numerical solution of dynamic portfolio optimization with transaction costs. Working paper. EID arXiv :2003.01809 . Chu, M. T., Kuo, C.-H., & Lin, M. M. (2013). Tensor spline approximation in economic dynamics with uncertainties. Computational Economics, 42(2), 175–198. Cocco, J. F., Gomes, F. J., & Maenhout, P. J. (2005). Consumption and portfolio choice over the life cycle. Review of Financial Studies, 18(2), 491–533. Constantinides, G. M. (1979). Multiperiod consumption and investment behavior with convex transactions costs. Management Science, 25(11), 1127–1137. Cox, M. G. (1972). The numerical evaluation of B-splines. IMA Journal of Applied Mathematics, 10(2), 134–149. de Boor, C. (1972). On calculating with B-splines. Journal of Approximation Theory, 6(1), 50–62. De Giorgi, E. G., & Legg, S. (2012). Dynamic portfolio choice and asset pricing with narrow framing and probability weighting. Journal of Economic Dynamics and Control, 36(7), 951–972. Epstein, L. G., & Zin, S. E. (1989). Substitution, risk aversion, and the temporal behavior of consumption and asset returns: A theoretical framework. Econometrica, 57(4), 937–969. Fletcher, R. (2013). Practical methods of optimization. New York: Wiley. https ://doi. org/10.1002/97811 18723 203. Garlappi, L., & Skoulakis, G. (2009). Numerical solutions to dynamic portfolio problems: The case for value function iteration using Taylor approximation. Computational Economics, 33(2), 193–207. Gill, P. E., Murray, W., & Saunders, M. A. (2005). SNOPT: An SQP algorithm for large-scale constrained optimization. SIAM Review, 47(1), 99–131. Habermann, C., & Kindermann, F. (2007). Multidimensional spline interpolation: Theory and applications. Computational Economics, 30(2), 153–169. Horneff, V., Maurer, R., & Schober, P. (2016). Efficient parallel solution methods for dynamic portfolio choice models in discrete time. Working paper 2665031. Available at SSRN. Goethe University Frankfurt. https ://doi.org/10.2139/ssrn.26650 31. 223 1 3 Solving High‑Dimensional Dynamic Portfolio Choice Models… Horneff, W. J., Maurer, R., & Rogalla, R. (2010). Dynamic portfolio choice with deferred annuities. Journal of Banking & Finance, 34(11), 2652–2664. Horneff, W. J., Maurer, R., & Stamos, M. Z. (2008). Life-cycle asset allocation with annuity markets. Journal of Economic Dynamics and Control, 32(11), 3590–3612. Hubener, A., Maurer, R., & Mitchell, O. S. (2016). How family status and social security claiming options shape optimal life cycle portfolios. Review of Financial Studies, 29(4), 937–978. Hubener, A., Maurer, R., & Rogalla, R. (2014). Optimal portfolio choice with annuities and life insurance for retired couples. Review of Finance, 18(1), 147–188. Höllig, K., & Hörner, J. (2013). Approximation and modeling with B-splines. SIAM. http://books tore. siam.org/ot132 /. Inkmann, J., Lopes, P., & Michaelides, A. (2011). How deep is the annuity market participation puzzle? Review of Financial Studies, 24(1), 279–319. Judd, K. L. (1998). Numerical methods in economics. Cambridge: MIT Press. https ://mitpr ess.mit. edu/books /numer icalm ethod s-econo mics. Judd, K. L., Maliar, L., Maliar, S., & Valero, R. (2014). Smolyak method for solving dynamic economic models: Lagrange interpolation, anisotropic grid and adaptive domain. Journal of Economic Dynamics and Control, 44, 92–123. Judd, K. L., & Solnick, A. (1994). Numerical dynamic programming with shape-preserving splines. Unpublished manuscript, Hoover Institution. https ://web.stanf ord.edu/~judd/paper s/dpsha pe.pdf. Kamin, J. H. (1975). Optimal portfolio revision with a proportional transaction cost. Management Science, 21(11), 1263–1271. Liu, H. (2004). Optimal consumption and investment with transaction costs and multiple risky assets. Journal of Finance, 59(1), 289–338. Liu, H., & Loewenstein, M. (2002). Optimal portfolio selection with transaction costs and finite horizons. Review of Financial Studies, 15(3), 805–835. Magill, M. J., & Constantinides, G. M. (1976). Portfolio selection with transactions costs. Journal of Economic Theory, 13(2), 245–263. Merton, R. C. (1969). Lifetime portfolio selection under uncertainty: The continuous-time case. Review of Economics and Statistics, 51(3), 247–257. Muthuraman, K., & Kumar, S. (2006). Multidimensional portfolio optimization with proportional transaction costs. Mathematical Finance, 16(2), 301–335. Pflüger, D. (2010). Spatially adaptive sparse grids for high-dimensional problems. Verlag Dr. Hut. https :// www5.in.tum.de/pub/pflue ger10 spati ally.pdf. Pflüger, D. (2012). Spatially adaptive refinement. In J. Garcke & M. Griebel (Eds.), Sparse grids and applications. Lecture Notes in Computational Science and Engineering (Vol. 88, pp. 243–262). Berlin: Springer. Philbrick, C. R, Jr., & Kitanidis, P. K. (2001). Improved dynamic programming methods for optimal control of lumped-parameter stochastic systems. Operations Research, 49(3), 398–412. Rust, J. (2008). Dynamic programming. In S. N. Durlauf & L. E. Blume (Eds.), The new Palgrave dictionary of economics (Vol. 1-8, pp. 1471–1489). London: Palgrave Macmillan. Schober, P. (2018). Solving dynamic portfolio choice models in discrete time using spatially adaptive sparse grids. In J. Garcke, D. Pflüger, C. Webster, & G. Zhang (Eds.), Sparse grids and applications-Miami 2016. Lecture Notes in Computational Science and Engineering (Vol. 123, pp. 135– 173). Berlin: Springer. Sickel, W., & Ullrich, T. (2011). Spline interpolation on sparse grids. Applicable Analysis, 90(3–4), 337–383. Stoyanov, M. (2017). User manual: TASMANIAN sparse grids v4.0. Technical report, Oak Ridge National Laboratory. https ://tasma nian.ornl.gov/docum ents/UserM anual .pdf. Valentin, J. (2019). B-splines for sparse grids: Algorithms and application to higher-dimensional optimization. Ph.D. thesis. University of Stuttgart. Valentin, J., & Pflüger, D. (2016). Hierarchical gradient-based optimization with B-splines on sparse grids. In J. Garcke & D. Pflüger (Eds.), Sparse grids and applications-Stuttgart 2014. Lecture Notes in Computational Science and Engineering (Vol. 109, pp. 315–336). Berlin: Springer. Valentin, J., Sprenger, M., Pflüger, D., & Röhrle, O. (2018). Gradient-based optimization with B-splines on sparse grids for solving forward-dynamics simulations of three-dimensional, continuum-mechanical musculoskeletal system models. International Journal for Numerical Methods in Biomedical Engineering, 34(5), 1–21. 224 P.Schober et al. 1 3 Winschel, V., & Krätzig, M. (2010). Solving, estimating, and selecting nonlinear dynamic models without the curse of dimensionality. Econometrica, 78(2), 803–821. Zenger, C. (1991). Sparse grids. In W. Hackbusch (Ed.), Parallel algorithms for partial differential equations. Notes on Numerical Fluid Mechanics (Vol. 31, pp. 241-251). Braunschweig: Vieweg. http:// www5.in.tum.de/pub/zenge r91sg .pdf. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.