$\displaystyle=\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}\left(-\frac{1}{\beta}\operatorname{Tr}{\left[\log{\left(1+e^{-\beta
H_{\text{BdG}}^{({\bm{\lambda}})}}\right)}\right]}\right),$ (216)
$\displaystyle=2\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\alpha}}\partial_{\lambda_{\beta}}\partial_{\lambda_{\gamma}}F_{\bm{\lambda}}.$
(217)

As a result, we reproduce the formula for the nonreciprocal superfluid density
effect given by

$2\sigma^{\alpha;\beta\gamma}_{\text{NRSF}}\left(\omega_{1}+\omega_{2};\omega_{1},\omega_{2}\right)=\frac{1}{2\omega_{1}\omega_{2}}A_{\alpha\beta\gamma}=\frac{1}{\omega_{1}\omega_{2}}f^{\alpha\beta\gamma}.$
(218)

The obtained formula is the generalized expression of Eq. (167) and captures
the nonreciprocal superfluid density effect in general NRO responses including
those other than the photocurrent generation $\omega_{1}=-\omega_{2}$.

Next, we consider the coefficients of the $O(\omega^{-1})$ term in Eq. (195).
Similarly to $A_{\alpha\beta\gamma}$, the coefficient $B_{\alpha\beta\gamma}$
is simplified as

$B_{\alpha\beta\gamma}=\int\frac{d\Omega}{2\pi
i}f(\Omega)\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\operatorname{Tr}{\left[{J}^{\alpha}_{{\bm{\lambda}}}\left(G_{\bm{\lambda}}^{R}(\Omega)-G_{\bm{\lambda}}^{A}(\Omega)\right){J}^{\gamma}_{{\bm{\lambda}}}F_{\bm{\lambda}}^{A}-{J}^{\alpha}_{{\bm{\lambda}}}F_{\bm{\lambda}}^{R}{J}^{\gamma}_{{\bm{\lambda}}}\left(G_{\bm{\lambda}}^{R}(\Omega)-G_{\bm{\lambda}}^{A}(\Omega)\right)\right]}.$
(219)

Since the $\lambda_{\beta}$-derivative is taken over the whole expression,
this contribution vanishes in the normal state. Interestingly, this formula
has the similar form to the regular part of the linear static conductivity.
The expression is given by

$\sigma_{\alpha\beta}^{({\bm{\lambda}})}=-i\int\frac{d\Omega}{2\pi
i}f(\Omega)\operatorname{Tr}{\left[{J}^{\alpha}_{{\bm{\lambda}}}F_{{\bm{\lambda}}}^{R}(\Omega){J}^{\beta}_{{\bm{\lambda}}}\left(G_{{\bm{\lambda}}}^{R}(\Omega)-G_{{\bm{\lambda}}}^{A}(\Omega)\right)-{J}^{\alpha}_{{\bm{\lambda}}}\left(G_{{\bm{\lambda}}}^{R}(\Omega)-G_{{\bm{\lambda}}}^{A}(\Omega)\right){J}^{\beta}_{{\bm{\lambda}}}F_{{\bm{\lambda}}}^{A}(\Omega)\right]}.$
(220)

Thus, Eq. (219) is rewritten as

$B_{\alpha\beta\gamma}=-i\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\sigma_{\alpha\gamma}^{({\bm{\lambda}})}.$
(221)

The superscript in $\sigma^{({\bm{\lambda}})}$ indicates that the Green
functions and current operators contain the parameter ${\bm{\lambda}}$. We
also recast $B_{\alpha\beta\gamma}^{\prime}$ as

$B^{\prime}_{\alpha\beta\gamma}=-i\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\gamma}}\sigma_{\alpha\beta}^{({\bm{\lambda}})}.$
(222)

We therefore call the low-frequency divergence from $B_{\alpha\beta\gamma}$
and $B^{\prime}_{\alpha\beta\gamma}$ (static) conductivity derivative effect.
The formula is given by

$2\sigma^{\alpha;\beta\gamma}_{\text{sCD}}=\frac{1}{2\omega_{2}}B_{\alpha\beta\gamma}+\frac{1}{2\omega_{1}}{B^{\prime}}_{\alpha\beta\gamma}=-\frac{i}{2}\lim_{\bm{\lambda}\rightarrow\bm{0}}\left(\frac{1}{\omega_{2}}\partial_{\lambda_{\beta}}\sigma_{\alpha\gamma}^{({\bm{\lambda}})}+\frac{1}{\omega_{1}}\partial_{\lambda_{\gamma}}\sigma_{\alpha\beta}^{({\bm{\lambda}})}\right).$
(223)

Finally, we present the symmetry analysis of the anomalous NRO conductivity in
Eq. (195). First, we discuss the nonreciprocal superfluid density effect in
the Green function representation. Taking Eq. (196), we apply the time-
reversal operation $\theta$ to what is in the parenthesis. Owing to the
antiunitary property of $\theta$, the expression is transformed as

$\displaystyle\operatorname{Tr}\left[{J}^{\alpha\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right]$
$\displaystyle=\operatorname{Tr}\left[\theta\left(G^{R}(\Omega)-G^{A}(\Omega)\right)^{\dagger}\theta^{-1}\theta{{J}^{\alpha\beta\gamma}}^{\dagger}\theta^{-1}\right],$
(224)
$\displaystyle=\operatorname{Tr}\left[\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\left(-{J}^{\alpha\beta\gamma}\right)\right],$
(225)
$\displaystyle=-\operatorname{Tr}\left[{J}^{\alpha\beta\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right],$
(226) $\displaystyle=0.$ (227)

The time-reversal operation is performed in the first line and the
$\mathcal{T}$-odd property of ${J}^{\alpha\beta\gamma}$ is implemented in the
second line. As a result, the nonreciprocal superfluid density is forbidden by
the $\mathcal{T}$ symmetry. It is similarly forbidden by the $\mathcal{P}$
symmetry, while the $\mathcal{PT}$ symmetry does not make any symmetry
constraint.

Following parallel calculations, we can classify $B_{\alpha\beta\gamma}$ by
the $\mathcal{T}$ and $\mathcal{PT}$ symmetries. Imposing the $\mathcal{T}$
symmetry, the first term in Eq. (202) is transformed as

$\operatorname{Tr}\left[{J}^{\alpha\beta}F^{R}(\Omega){J}^{\gamma}\left(G^{R}(\Omega)-G^{A}(\Omega)\right)\right]=-\operatorname{Tr}\left[{J}^{\alpha\beta}\left(G^{R}(\Omega)-G^{A}(\Omega)\right){J}^{\gamma}F^{R}\right],$
(228)

since the $\mathcal{T}$-parities of the paramagnetic and diamagnetic current
operators are odd and even, respectively. Performing similar symmetry
operations for other terms of $B_{\alpha\beta\gamma}$, we obtain the
$\mathcal{T}$ symmetric component given by the antisymetric part of the linear
conductivity tensor

$B_{\alpha\beta\gamma}^{(\mathcal{T})}=-\frac{i}{2}\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\left(\sigma_{\alpha\gamma}^{({\bm{\lambda}})}-\sigma_{\gamma\alpha}^{({\bm{\lambda}})}\right).$
(229)

Since the
$\sigma_{\alpha\gamma}^{({\bm{\lambda}})}-\sigma_{\gamma\alpha}^{({\bm{\lambda}})}$
is expressed by the Berry curvature, this anomalous NRO response is given by
the Berry curvature derivative. We reproduce the expression for the Berry
curvature derivative by performing the trace. Using the eigenstates of the BdG
Hamiltonian $H_{\text{BdG}}^{({\bm{\lambda}})}$, we get the relation

$\displaystyle\int\frac{d\Omega}{2\pi
i}f(\Omega)\operatorname{Tr}{\left[-{J}_{\bm{\lambda}}^{\alpha}F_{{\bm{\lambda}}}^{R}(\Omega){J}_{\bm{\lambda}}^{\gamma}\left(G_{\bm{\lambda}}^{R}(\Omega)-G_{\bm{\lambda}}^{A}(\Omega)\right)+{J}_{\bm{\lambda}}^{\alpha}\left(G_{\bm{\lambda}}^{R}(\Omega)-G_{\bm{\lambda}}^{A}(\Omega)\right){J}_{\bm{\lambda}}^{\gamma}F_{\bm{\lambda}}^{R}(\Omega)\right)}$
$\displaystyle=\sum_{a\neq
b}\frac{f_{ab}}{E_{ab}^{2}}\left({J}_{\bm{\lambda}}^{\alpha}\right)_{ab}\left({J}_{\bm{\lambda}}^{\gamma}\right)_{ba},$
(230) $\displaystyle=\sum_{a\neq
b}f_{a}\left[\left(\xi^{\lambda_{\alpha}}\right)_{ab}\left(\xi^{\lambda_{\gamma}}\right)_{ba}-\left(\alpha\leftrightarrow\gamma\right)\right],$
(231)

and Eq. (231) is recast as

$\sum_{a\neq
b}f_{a}\left[\left(\xi^{\lambda_{\alpha}}\right)_{ab}\left(\xi^{\lambda_{\gamma}}\right)_{ba}-\left(\alpha\leftrightarrow\gamma\right)\right]=\sum_{a}-i\epsilon_{\alpha\gamma\beta}\Omega_{a}^{\beta}f_{a}.$
(232)

Finally, we rewrite the static conductivity derivative effect for
$\mathcal{T}$ symmetric case by

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{sCD$(\mathcal{T})$}}$
$\displaystyle=\frac{1}{2\omega_{2}}B^{(\mathcal{T})}_{\alpha\beta\gamma}+\frac{1}{2\omega_{1}}{B^{\prime}}^{(\mathcal{T})}_{\alpha\beta\gamma},$
(233)
$\displaystyle=-\frac{i}{2}\lim_{\bm{\lambda}\rightarrow\bm{0}}\left[\frac{1}{\omega_{2}}\epsilon_{\alpha\gamma\delta}\partial_{\lambda_{\beta}}\left(\sum_{a}\Omega^{\delta}_{a}f_{a}\right)+\frac{1}{\omega_{1}}\epsilon_{\alpha\beta\delta}\partial_{\lambda_{\gamma}}\left(\sum_{a}\Omega^{\delta}_{a}f_{a}\right)\right],$
(234)
$\displaystyle=-\frac{i}{2}\left(\epsilon_{\alpha\beta\delta}\frac{B_{d}^{\gamma\delta}}{\omega_{1}}+\epsilon_{\alpha\gamma\delta}\frac{B_{d}^{\beta\delta}}{\omega_{2}}\right),$
(235)

in the low-frequency regime.

Next, we consider the symmetry constraint of the $\mathcal{PT}$ symmetry on
the coefficient $B_{\alpha\beta\gamma}$. On the basis of the parallel
discussion, the antisymmetric component is forbidden. Accordingly,
$B_{\alpha\beta\gamma}$ is expressed by the symmetric component

$B_{\alpha\beta\gamma}^{(\mathcal{PT})}=-\frac{i}{2}\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\left(\sigma_{\alpha\gamma}^{({\bm{\lambda}})}+\sigma_{\gamma\alpha}^{({\bm{\lambda}})}\right),$
(236)

which leads to the Drude derivative $\hat{D_{d}}$ defined by

$D_{d}^{\beta;\alpha\gamma}=iB_{\alpha\beta\gamma}^{(\mathcal{PT})}=\lim_{\bm{\lambda}\rightarrow\bm{0}}\partial_{\lambda_{\beta}}\left(\int\frac{d\Omega}{2\pi}\left(-\frac{\partial
f(\Omega)}{\partial\Omega}\right)\sum_{a}\left({J}_{\bm{\lambda}}^{\alpha}\right)_{aa}\left({J}_{\bm{\lambda}}^{\gamma}\right)_{aa}2\Im{\left[\left(G_{\bm{\lambda}}^{R}(\Omega)\right)_{aa}\right]}^{2}\right).$
(237)

This represents the intraband effect. The Drude derivative $\hat{D}_{d}$
vanishes in the absence of the Bogoliubov quasiparticles since it contains the
Fermi-surface factor $\partial_{\Omega}f(\Omega)$. Note that the parameter
${\bm{\lambda}}$ does not break the $\mathcal{PT}$ symmetry and we can
directly apply the $\mathcal{PT}$ symmetry to Eq. (219). The obtained result
is the same as Eq. (237). As a result, we obtain the static conductivity
derivative effect for $\mathcal{PT}$ symmetric case as

$\displaystyle 2\sigma^{\alpha;\beta\gamma}_{\text{sCD$(\mathcal{PT})$}}$
$\displaystyle=\frac{1}{2\omega_{2}}B^{(\mathcal{PT})}_{\alpha\beta\gamma}+\frac{1}{2\omega_{1}}{B^{\prime}}^{(\mathcal{PT})}_{\alpha\beta\gamma},$
(238)
$\displaystyle=-\frac{i}{2\omega_{2}}D_{d}^{\beta;\alpha\gamma}-\frac{i}{2\omega_{1}}D_{d}^{\gamma;\alpha\beta}.$
(239)

## References

  * (1) M. Fiebig, V. V. Pavlov, and R. V. Pisarev, J. Opt. Soc. Am. B 22, 96 (2005).
  * (2) J. Orenstein, J. E. Moore, T. Morimoto, D. H. Torchinsky, J. W. Harter, and D. Hsieh, Annu. Rev. Condens. Matter Phys. 12, 247 (2021).
  * (3) L. Zhao, C. A. Belvin, R. Liang, D. A. Bonn, W. N. Hardy, N. P. Armitage, and D. Hsieh, Nat. Phys. 13, 250 (2017).
  * (4) A. de la Torre, K. L. Seyler, L. Zhao, S. Di Matteo, M. S. Scheurer, Y. Li, B. Yu, M. Greven, S. Sachdev, M. R. Norman, and D. Hsieh, arXiv:2008.06516.
  * (5) B. I. Sturman and V. M. Fridkin, Photovoltaic and Photorefractive Effects in Noncentrosymmetric Materials (Gordon and Breach, New York, 1992).
  * (6) J. Liu, F. Xia, D. Xiao, F. J. García de Abajo, and D. Sun, Nat. Mater. 19, 830 (2020).
  * (7) N. Ogawa, R. Yoshimi, K. Yasuda, A. Tsukazaki, M. Kawasaki, and Y. Tokura, Nat. Commun. 7, 12246 (2016).
  * (8) A. M. Burger, L. Gao, R. Agarwal, A. Aprelev, J. E. Spanier, A. M. Rappe, and V. M. Fridkin, Phys. Rev. B 102, 081113 (2020).
  * (9) M. Tinkham, Introduction to superconductivity (McGraw-Hill, New York, 1996).
  * (10) Y. Tokura and N. Nagaosa, Nat. Commun. 9, 3740 (2018).
  * (11) T. Ideue, S. Koshikawa, H. Namiki, T. Sasagawa, and Y. Iwasa, Phys. Rev. Research 2, 042046 (2020).
  * (12) F. Ando, Y. Miyasaka, T. Li, J. Ishizuka, T. Arakawa, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Nature 584, 373 (2020).
  * (13) T. Xu, T. Morimoto, and J. E. Moore, Phys. Rev. B 100, 220501 (2019).
  * (14) X. Yang, C. Vaswani, C. Sundahl, M. Mootz, L. Luo, J. H. Kang, I. E. Perakis, C. B. Eom, and J. Wang, Nat. Photonics 13, 707 (2019).
  * (15) S. Nakamura, K. Katsumi, H. Terai, and R. Shimano, Phys. Rev. Lett. 125, 097004 (2020).
  * (16) C. Vaswani, M. Mootz, C. Sundahl, D. H. Mudiyanselage, J. H. Kang, X. Yang, D. Cheng, C. Huang, R. H. J. Kim, Z. Liu, L. Luo, I. E. Perakis, C. B. Eom, and J. Wang, Phys. Rev. Lett. 124, 207003 (2020).
  * (17) S. Lim, C. M. Varma, H. Eisaki, and A. Kapitulnik, arXiv:2011.06755.
  * (18) D. Aoki, K. Ishida, and J. Flouquet, J. Phys. Soc. Japan 88, 022001 (2019).
  * (19) J. A. Sauls, Adv. Phys. 43, 113 (1994).
  * (20) R. Joynt and L. Taillefer, Rev. Mod. Phys. 74, 235 (2002).
  * (21) H. Watanabe and Y. Yanase, Phys. Rev. X 11, 011001 (2021).
  * (22) E. Bauer and M. Sigrist (eds.), Non-Centrosymmetric Superconductors: Introduction and Overview, Vol. 847 (Springer, Berlin, 2012).
  * (23) Y. Saito, Y. Nakamura, M. S. Bahramy, Y. Kohama, J. Ye, Y. Kasahara, Y. Nakagawa, M. Onga, M. Tokunaga, T. Nojima, Y. Yanase, and Y. Iwasa, Nat. Phys. 12, 144 (2016).
  * (24) J. M. Lu, O. Zheliuk, I. Leermakers, N. F. Q. Yuan, U. Zeitler, K. T. Law, and J. T. Ye, Science 350, 1353 (2015).
  * (25) T. Schumann, L. Galletti, H. Jeong, K. Ahadi, W. M. Strickland, S. Salmani-Rezaie, and S. Stemmer, Phys. Rev. B 101, 100503(R) (2020).
  * (26) J. T. Ye, Y. J. Zhang, R. Akashi, M. S. Bahramy, R. Arita, and Y. Iwasa, Science 338, 1193 (2012).
  * (27) N. Reyren, S. Thiel, A. D. Caviglia, L. F. Kourkoutis, G. Hammerl, C. Richter, C. W. Schneider, T. Kopp, A.-S. Rüetschi, D. Jaccard, M. Gabay, D. A. Muller, J.-M. Triscone, and J. Mannhart, Science 317, 1196 (2007).
  * (28) Z. Hiroi, J.-I. Yamaura, T. C. Kobayashi, Y. Matsubayashi, and D. Hirai, J. Phys. Soc. Japan 87, 024702 (2018).
  * (29) S. Sumita, T. Nomoto, and Y. Yanase, Phys. Rev. Lett. 119, 027001 (2017).
  * (30) P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
  * (31) Y. Wang and L. Fu, Phys. Rev. Lett. 119, 187003 (2017).
  * (32) S. Ran, C. Eckberg, Q. P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I. L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Science 365, 684 (2019).
  * (33) S. Ran, I.-L. Liu, Y. S. Eo, D. J. Campbell, P. M. Neves, W. T. Fuhrman, S. R. Saha, C. Eckberg, H. Kim, D. Graf, F. Balakirev, J. Singleton, J. Paglione, and N. P. Butch, Nat. Phys. 15, 1250 (2019).
  * (34) D. Braithwaite, M. Vališka, G. Knebel, G. Lapertot, J. P. Brison, A. Pourret, M. E. Zhitomirsky, J. Flouquet, F. Honda, and D. Aoki, Commun. Phys. 2, 1 (2019).
  * (35) D. Aoki, F. Honda, G. Knebel, D. Braithwaite, A. Nakamura, D. Li, Y. Homma, Y. Shimizu, Y. J. Sato, J.-P. Brison, and J. Flouquet, J. Phys. Soc. Japan 89, 053705 (2020).
  * (36) G. Knebel, M. Kimata, M. Vališka, F. Honda, D. Li, D. Braithwaite, G. Lapertot, W. Knafo, A. Pourret, Y. J. Sato, Y. Shimizu, T. Kihara, J.-P. Brison, J. Flouquet, and D. Aoki, J. Phys. Soc. Japan 89, 053707 (2020).
  * (37) W. C. Lin, D. J. Campbell, S. Ran, I. L. Liu, H. Kim, A. H. Nevidomskyy, D. Graf, N. P. Butch, and J. Paglione, npj Quantum Materials 5, 1 (2020).
  * (38) D. Aoki, M. Kimata, Y. J. Sato, G. Knebel, F. Honda, A. Nakamura, D. Li, Y. Homma, Y. Shimizu, W. Knafo, D. Braithwaite, M. Valiska, A. Pourret, J.-P. Brison, and J. Flouquet, arXiv:2104.01565.
  * (39) J. Ishizuka and Y. Yanase, Phys. Rev. B 103, 094504 (2021).
  * (40) C. Kallin and J. Berlinsky, Rep. Prog. Phys. 79, 054502 (2016).
  * (41) P. M. R. Brydon, D. S. L. Abergel, D. F. Agterberg, and V. M. Yakovenko, Phys. Rev. X 9, 031025 (2019).
  * (42) S. Yonezawa, Condensed Matter 4, 1 (2019).
  * (43) Y. Michishita and R. Peters, Phys. Rev. B 103, 195133 (2021).
  * (44) D. J. Passos, G. B. Ventura, J. M. V. P. Lopes, J. M. B. L. D. Santos, and N. M. R. Peres, Phys. Rev. B 97, 235446(2018).
  * (45) D. Vanderbilt, Berry Phases in Electronic Structure Theory: Electric Polarization,Orbital Magnetization and Topological Insulators (Cambridge University Press, Cambridge, 2018).
  * (46) R. von Baltz and W. Kraut, Phys. Rev. B 23, 5590 (1981).
  * (47) J. E. Sipe and A. I. Shkrebtii, Phys. Rev. B 61, 5337 (2000).
  * (48) F. de Juan, Y. Zhang, T. Morimoto, Y. Sun, J. E. Moore, and A. G. Grushin, Phys. Rev. Research 2, 012017 (2020).
  * (49) W. Kohn, Phys. Rev. 133, A171 (1964).
  * (50) R. Resta, J. Phys. Condens. Matter 30, 414001 (2018).
  * (51) Y. Zhang, T. Holder, H. Ishizuka, F. de Juan, N. Nagaosa, C. Felser, and B. Yan, Nat. Commun. 10, 3783 (2019).
  * (52) T. Holder, D. Kaplan, and B. Yan, Phys. Rev. Research 2, 033100 (2020).
  * (53) Q. Niu and D. J. Thouless, J. Phys. A 17, 2453 (1984).
  * (54) P. S. Halasyamani and K. R. Poeppelmeier, Chem. Mater. 10, 2753 (1998).
  * (55) J. E. Moore and J. Orenstein, Phys. Rev. Lett. 105, 026805 (2010).
  * (56) J. Ahn, G.-Y. Guo, and N. Nagaosa, Phys. Rev. X 10, 041041 (2020).
  * (57) D. J. Thouless, Phys. Rev. B 27, 6083 (1983).
  * (58) H. Watanabe and M. Oshikawa, Phys. Rev. B 102, 165137 (2020); H. Watanabe, Y. Liu, and M. Oshikawa, J. Stat. Phys. 181, 2050 (2020).
  * (59) D. J. Scalapino, S. R. White, and S-C. Zhang, Phys. Rev. B 47, 7995 (1993).
  * (60) Similar expression corresponding to Eq. (48) has recently been obtained in Ref. Resta2021  where the normal state is assumed and thus the obtained term may vanish.
  * (61) R. Resta, arXiv:2101.10949.
  * (62) R. Shimano and N. Tsuji, Annu. Rev. Condens. Matter Phys. 11, 103 (2020), and references therein.
  * (63) H. Watanabe, A. Daido, and Y. Yanase, submitted.
  * (64) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * (65) Y. Yanase, J. Phys. Soc. Japan 83, 014703 (2014).
  * (66) J. Železný, H. Gao, K. Výborný, J. Zemen, J. Mašek, A. Manchon, J. Wunderlich, J. Sinova, and T. Jungwirth, Phys. Rev. Lett. 113, 157201 (2014).
  * (67) S. Harada, J. J. Zhou, Y. G. Yao, Y. Inada, and G.-q. Zheng, Phys. Rev. B 86, 220502 (2012).
  * (68) S. Kitagawa, K. Ishida, T. C. Kobayashi, Y. Matsubayashi, D. Hirai, and Z. Hiroi, J. Phys. Soc. Japan 89, 053701 (2020).
  * (69) L. Gao, Z. Addison, E. J. Mele, and A. M. Rappe, arXiv.2011.06542.
  * (70) Y. Gao and D. Xiao, Phys. Rev. Lett. 122, 227402 (2019).
  * (71) K. Nomura, S. Ryu, A. Furusaki, and N. Nagaosa, Phys. Rev. Lett. 108, 026802 (2012).
  * (72) H. Sumiyoshi and S. Fujimoto, J. Phys. Soc. Japan 82, 023602 (2013).
  * (73) M. Tinkham and R. A. Ferrell, Phys. Rev. Lett. 2, 331 (1959).
  * (74) F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore, Nat. Commun. 8, 15995 (2017).
  * (75) E. Deyo, L. E. Golub, E. L. Ivchenko, and B. Spivak, arXiv.0904.1917.
  * (76) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * (77) D. E. Parker, T. Morimoto, J. Orenstein, and J. E. Moore, Phys. Rev. B 99, 045121 (2019).


# Optimal liquidation with temporary and permanent price impact, an
application to cryptocurrencies

Hugo E. Ramirez<EMAIL_ADDRESS>Universidad del Rosario, Calle
12C No. 4-69, Bogotá, Colombia  Julián Fernando Sánchez López
<EMAIL_ADDRESS>Universidad del Rosario, Calle 12C No. 4-69,
Bogotá, Colombia

###### Abstract

This paper studies the optimal liquidation of stocks in the presence of
temporary and permanent price impacts, and we focus in the case of
cryptocurrencies. We start by presenting analytical solutions to the problem
with linear temporary impact, and linear and quadratic permanent impact. Then,
using data from the order book of the BNB cryptocurrency, we estimate the
functional form of the temporary and permanent price impact in three different
scenarios: underestimation, overestimation and average estimation, finding
different functional forms for each scenario. Using finite differences and
optimal policy iteration, we solve the problem numerically and observe
interesting changes in the optimal liquidation policy when applying calibrated
linear and power forms for the temporary and permanent price impacts. Then,
with these optimal policies, we identify optimal liquidation trajectories and
simulate the liquidation of initial inventories to compare the performance
among the optimal strategies under different parametrizations and against a
naive strategy. Finally, we characterize the optimal policies based on the
functional form of the inventory and find that policies generating the highest
revenue are those starting with a low trading rate and increasing it as time
passes.

Keywords— Optimal liquidation, price impact, finite differences, Stochastic
Control, cryptocurrencies

##  1 Introduction

The market of cryptocurrencies is increasing in importance since the net value
of the market has grown from $\$3.8$ billion $USD$ in $2015$ to $\$943.41$
billion $USD$ in the second quarter of $2022$, having a maximum of over
$\$3000$ Billion $USD$111Data from
https://www.statista.com/statistics/730876/cryptocurrency-maket-value/ and
https://coinmarketcap.com/. Likewise, the use of cryptocurrencies has
advantages such as its articulation with Decentralized Finance (DeFi)222In
accordance with Piñeiro-Chousa et al.,  (2022) DeFi can be understood as the
process that makes use of the blockchain for the development and
implementation of novel financial products and services., security and
transparency, among others, which, according to Blanco Encinosa,  (2021), make
them very convenient investment alternatives. However, this type of investment
is very volatile and many investors may have or want to acquire or liquidate
their positions optimally on a given time period. In this paper, we focus on
optimal strategies for liquidation of assets considering price impacts,
specifically different forms of temporary and permanent price impact
calibrated from cryptocurrencies.

Optimal liquidation consists in finding the optimal strategy an agent must
follow when selling large amounts of shares over a finite period of time and
minimizing adverse effects, consequences of her own actions. The time
dependence is of great importance because, if the agent executes in a short
time, by selling fast, the rate of execution increases and so does the impact
on the price of shares. But if she sells at a slow pace over a long period of
time, holdings are exposed to a greater uncertainty due to the volatility,
which is crucial in the case of cryptocurrencies. Thus, in our setting, the
agent must sell a number $\Upsilon$ of shares in a given time frame $[0,T)$,
where she must find the optimal selling rate $\nu_{t}$ such that the revenue
from liquidation is maximized. To this end, we adopt an optimization problem
based on a model that includes linear and non-linear impact functions, and is
calibrated with data from the Limit Order Book (LOB) of the BNB
cryptocurrency. To tackle the problem, we use stochastic control techniques
and present two original analytical solutions in specific scenarios, but in
more general settings and when analytical solutions were not at hand, we use
numerical methods, that is finite differences over a non-linear Partial
Differential Equation (PDE).

The well-known model for optimal execution of transactions in Almgren and
Chriss,  (2000) establishes two types of impact for the asset price: temporary
and permanent. In their work, temporary impact refers to “…temporary
imbalances in supply in demand caused by our trading leading to temporary
price movements away from equilibrium” and permanent impact “…means changes in
the equilibrium price due to our trading which remain at least for the life of
our liquidation”. Theoretical approaches agree that the permanent price impact
must be linear, Gatheral,  (2010) shows this condition is necessary to avoid
dynamic arbitrage. He presents the principle of no-dynamic-arbitrage which
states that the trading cost is non-negative for any round-trip trade
strategy, i.e. $\displaystyle\int_{0}^{T}\nu_{t}dt=0$, and shows that non-
linear permanent market impact is inconsistent with the principle of no-
dynamic-arbitrage. Conversely, our study focuses in only liquidating shares,
and then the sum of this sequence of trades cannot be zero, therefore the
principle of no-dynamic arbitrage is not applicable. Consequently our model
considers non-linear permanent price impacts, which are estimated from data,
and studies the effect of such impact functions on the optimized trading rate.

Additionally, empirical evidence shows that market impact is not linear, but
similar to a square root. That is, it is proportional to the square root of
the volume of shares executed, as seen in Toth et al.,  (2016) and Almgren et
al.,  (2005), or to the square root of the trade duration, as in Bershova and
Rakhlin,  (2013). Moreover, different approaches using nonlinear temporary or
permanent price impact have been studied in the literature, for example
Gueant,  (2014) relates to the result in Gatheral,  (2010) and extends some
theoretical results to the nonlinear setting, Barger and Lorig,  (2018) offers
a theoretical approach, considering stochastic temporary and permanent price
impacts. Alfonsi and Schied,  (2010) proposes a model where the impact depends
on the theoretical shape of LOB given by a density (exponential), whereas our
model does not rely on any shape and is thus capable of capture the actual
book in simulations. More recently Brunovský et al.,  (2018) solves for
optimal liquidation with linear price impact, although importantly shows some
theoretical advances by adding the possibility of acquiring shares as well as
selling.

Although related literature shows different perspectives to the solution of
the optimal liquidation problem for large amounts of shares, in this paper we
study the problem based on the setting proposed in Cartea et al.,  (2015), but
using more general functional forms for temporary and permanent price impact
(TPI and PPI respectively). At first, we extend the theoretical model in
Cartea et al.,  (2015) by considering linear TPI, and linear and quadratic
PPI, and show close form solutions for these models. Since analytical
solutions are scarce, for more general models we use numerical methods,
specifically finite differences and optimal policy iteration, and solve the
problem for different forms of the TPI and PPI functions, i.e. functions
result of calibration. This type of problems, i.e. a PDE with two dimensions
in space and only one diffusion, are frequently numerically solved using the
ADI scheme, but have some stability issues, as shown in the classical work of
Peaceman and Ratchford Jr,  (1955), henceforth also, for example, in Thomas,
(1995),Tavella and Randall,  (2000) and more recently by Duffy,  (2006). So in
order to improve convergence in the multi-dimensional finite differences we
propose the use of an alternative scheme, namely implicit directional scheme.

After that, using the technique in Cartea and Jaimungal, 2016a we calibrate,
in a cryptocurrency market, the functional forms of the TPI and PPI functions.
A first aspect to highlight is that calibrated parameters of the functional
forms for the TPI and PPI change according to the depth in ticks and therefore
the volumes of the LOB 333We refer to the volumes of the LOB as the number of
shares posted at each tick level. That is, given a fixed LOB setting (depth in
ticks and volumes), if the inventory to liquidate is sufficiently small, the
liquidation can be done barely walking the book, so that the impact would be
small and increasing as the liquidation target is reached, this configures a
power form $x^{a}$, where $a>1$. The opposite case occurs when the inventory
to be liquidated is large for the LOB setting, since at first the impact would
be large, but decreasing as the liquidation continues, this configures a power
form $x^{a}$ but having $0<a<1$. As a third scenario, we find a in-between
point, where the TPI and PPI approaches a linear form, i.e. $x^{a}$ with
$a=1$. Using this three configurations for calibration gives different
functions for TPI and PPI, we are able to identify the different optimal
trading policies for each of these scenarios. Thus, in accordance with Lillo
et al.,  (2003) and Almgren et al.,  (2005), modeling the problem must include
power forms for TPI and PPI, but for these more complex and general scenarios,
the use numerical techniques to find solutions is compulsory.

We find that optimal liquidation strategies change for different TPI and PPI
functions, and that these outperform a naive strategy of liquidating all the
inventory at a single point in time, which we confirm with simulations. We
observe that, we may describe these strategies in a common functional form of
the optimal inventory, as a non-increasing monotonic function $\displaystyle
q(t)=\frac{-\Upsilon}{T^{d_{2}}}t^{d_{2}}+\Upsilon$ subject to $q(0)=\Upsilon$
and $q(T)=0$. Under certain conditions and using real life data simulations,
similar to a paper trading, the optimal policies starting with a slow rate and
increasing it as time passes have the highest cumulative revenue. Thus, our
results in some sense confirm the conclusions about concavity of market impact
in Curato et al.,  (2017), stating that concave functions generate more
profit. This suggests that beyond the results in calibration, the performance
of optimal policies does not depend strongly on the functional form of TPI and
PPI, and that such performance may be characterized by the initial amount of
shares to be liquidated compared with the availability of the LOB. This has an
important potential in practice, since eliminates the hustle of calibrating
models and thus saves time by implementing directly the optimal strategy.

This paper is structured as follows: section $2$ introduces the model, based
on Cartea et al.,  (2015), and analytical solutions of the problem of
liquidation from the perspective of optimal stochastic control. Section $3$
presents analytical solutions when the TPI is linear, and PPI is linear and
quadratic as extensions of Cartea et al.,  (2015), and a general solution
using finite differences, for more general TPI and PPI. Section $4$ shows the
estimation of a possible functional forms of TPI and PPI and to calibration of
the parameters for the BNB cryptocurrency, using the LOB data. In section $5$,
we find the optimal policies and analyze them by varying some of the
calibrated parameters. Also, we compare the performance of numeric and naive
trading strategies simulating the optimal liquidation of initial inventories
on LOB data. Finally, we establish a characterization of the optimal policies
from the optimal inventory available, and independent of the calibration. In
section $6$ we give some final remarks.

##  2 The model

We work on a probability space $\left(\Omega,\mathcal{F},\mathbb{P}\right)$
and build on the model proposed by Cartea et al.,  (2015), the agent has an
amount $\Upsilon$ of shares of an asset that she wants to liquidate in the
time interval $\left[0,T\right)$ and obtain the most profit. The inventory
changes through time and has dynamics ${dQ}^{\nu}_{t}={-\nu}_{t}dt$, thus
$\nu_{t}$ is the trade rate. The mid price of the asset is governed by the SDE

$\displaystyle{dS}^{\nu}_{t}=-g\left(\nu_{t}\right)dt+\sigma dW_{t},$ (1)

where $W_{t}$ is a standard Brownian motion, $\sigma\in\mathbb{R}^{+}$
represents the volatility and $g:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is the
permanent price impact, dependent of the trade rate. When an agent sells
shares of the stock the actual price of the transaction is called the
execution price and follows

$\displaystyle{\hat{S}}^{\nu}_{t}=\left(S^{\nu}_{t}-\frac{1}{2}\Delta\right)-f\left(\nu_{t}\right),$
(2)

where $\Delta$ is the bid-ask spread in the LOB and the function
$f:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is the temporary price impact. Finally the
agent wants to maximize her final expected utility wealth, earnings from
transactions, that is

$\displaystyle\mathbb{E}\left\\{\int^{T}_{0}{{\hat{S}}^{\nu}_{t}\nu_{t}dt}\right\\}$
(3)

the rate $\nu_{t}$ is controlled by the agent at each time point $t$ and her
actions affect the asset’s liquidation price, thus, the optimization problem
can be represented by the agent’s value function

$H\left(t,S,q\right)=\sup_{\nu_{t}\in\mathcal{A}(t,T)}{\mathbb{E}\left\\{\int^{T}_{t}\left(\left(S^{\nu}_{r}-\frac{1}{2}\Delta\right)-f\left(\nu_{r}\right)\right)\nu_{r}dr\right\\}}$
(4)

where $\mathcal{A}(t,T)$ represents the admissible set of non-negative bounded
strategies and the equation (4) satisfies the
$\mathcal{H}\mathcal{J}\mathcal{B}$ partial differential equation (PDE) 444For
a detailed explanation, refer to chapter $6$ of Cartea et al.,  (2015)

$\displaystyle
0={\partial}_{t}H+\frac{1}{2}{\sigma}^{2}{\partial}_{ss}H+\sup_{\nu\in\mathcal{A}}\left\\{-{g(\nu)\partial}_{s}H-\nu{\partial}_{q}H+\left(S-\frac{1}{2}\Delta-f(\nu)\right)\nu\right\\}$
(5)

subject to

$\displaystyle\displaystyle H(T,S,q)$ $\displaystyle\rightarrow$
$\displaystyle-\infty,\text{ when }t\rightarrow T\text{ and }q>0,$ (6)
$\displaystyle\displaystyle H(t,S,0)$ $\displaystyle\rightarrow$
$\displaystyle 0,$ (7)

the first condition corresponds to the penalty, negative revenue, for
approaching T with positive inventory, and the second is complementary and
guarantees that at time $t$ the revenue is zero if the inventory is zero.

Next, we extend the model in Cartea et al.,  (2015) by changing the PPI to a
linear and quadratic function of the trade rate $\nu$. Although linear PPI is
treated in Cartea et al.,  (2015) and Cartea and Jaimungal, 2016b , our
approach is essentially different because of the type of penalty, in the range
$[0,T)$. Since our penalty is a boundary condition going to $-\infty$, we are
not allowing the manager to have any inventory at maturity. On the other hand,
Cartea et al.,  (2015) and Cartea and Jaimungal, 2016b have a quadratic
penalty, which, although costly, allows the manager to end up with some
inventory.

##  3 Model solution approaches

###  3.1 Analytical approach

It is always satisfactory to have analytical solutions, this section shows the
scenarios where we were able to extend a close form formula for solutions. The
following lemmas show the analytical solutions for the
$\mathcal{H}\mathcal{J}\mathcal{B}$ PDE (5) subject to conditions (6) and (7),
when PPI is linear and quadratic, respectively. These lemmas may be seen as
extensions of the solution given in section $6.3$ of Cartea et al.,  (2015).

######  Lemma 1 (Linear permanent price impact).

Let TPI and PPI be linear functions, that is
$f\left(\nu_{t}\right)=a_{1}\nu_{t}+a_{2}\ $, and $g(\nu_{t})=b_{1}\nu_{t}$
555For simplicity and to find an analytical solution of
$\mathcal{H}\mathcal{J}\mathcal{B}$ PDS (5) we use $g(\nu_{t})=b_{1}\nu_{t}$
instead of $g(\nu_{t})=b_{1}\nu_{t}+b_{2}$ since the independent term $b_{2}$
does not allow us to find an analytical solution. respectively, then the
solution for the $\mathcal{H}\mathcal{J}\mathcal{B}$ differential equation (5)
subject to (6) and (7) is

$H\left(t,S,q\right)=q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\left(\frac{b_{1}}{2}+\frac{a_{1}}{T-t}\right)q^{2},$ (8)

and the optimal strategy for the liquidation of $q_{t}$ shares in the time
interval $[t,T)$ is

$\nu_{t}^{*}=\frac{q_{t}}{T-t}.$ (9)

###### Proof.

See appendix A. ∎

######  Lemma 2 (Quadratic permanent price impact).

Let us have a linear TPI and a quadratic PPI, that is
$f\left(\nu_{t}\right)=a_{1}\nu_{t}+a_{2}$, and
$g(\nu_{t})=c_{1}\nu_{t}^{2}+c_{2}\nu_{t}\ $, respectively, then the
analytical solution for the $\mathcal{H}\mathcal{J}\mathcal{B}$ PDE (5)
subject to (6) and (7) is

$H\left(t,S,q\right)=\left\\{\begin{array}[]{lcc}q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\frac{c_{2}}{2}q^{2}-\frac{4\left(c_{1}q+a_{1}\right)^{3}}{9c_{1}^{2}(T-t)}&if&q\neq
0\\\ \\\ 0&if&q=0,\end{array}\right.$ (10)

and the optimal strategy for the liquidation of $q_{t}$ shares in the time
interval $[t,T)$ is

$\nu_{t}^{*}=\frac{2\left(c_{1}q_{t}+a_{1}\right)}{3c_{1}\left(T-t\right)}$
(11)

###### Proof.

See appendix B. ∎

######  Remark.

As in lemma 1 no independent term is included in the quadratic PPI function,
because this allows us calculate an analytical solution of the equation (5).

From lemmas 1 and 2 we appreciate that whereas $H$ is quadratic in the $q$
variable for the former, it is cubic for the latter, although linear for $S$
in both cases. This implies that the change of the PPI from linear to
quadratic affects in a similar manner the optimal revenue $H$. Although
interestingly note that the optimal policy in both cases is a linear function
of $\displaystyle\frac{q}{T-t}$.

Since obtaining analytical solutions to any type of PPI function is a hard
task, and in many cases even impossible, to extend our model, we use numerical
solutions. Specifically, we want to test optimal strategies in a
cryptocurrency market and use, besides linear, calibrated power TPI and PPI
functions for which we do not have analytical solutions.

###  3.2 Numerical approach

####  3.2.1 Finite differences

We use implicit finite differences to find the solution of the PDE (5). This
technique approximates the partial derivatives through discrete expressions
constructed from Taylor expansions and thereby giving an approximate solution
of the differential equation.

We denote the discretized version of the function $H$ by
$H_{i,j}^{k}=H\left(k\Delta t,i\Delta s,j\Delta q\right)$, where $t=k\Delta
t$, $s=i\Delta s$ and $q=j\Delta q$, additionally we set $S_{max}=N_{s}\Delta
s$, $Q_{max}=N_{q}\Delta q$ and $T=N_{t}\Delta t$. As boundary conditions are
great of importance in this technique, we state these next:

  * •

When the price is $S=0$, we think of it as an absorbent state and the possible
revenue for potential trades at this price is $0$, thus

$\displaystyle H(t,0,q)=0.$ (12)

  * •

If the inventory reaches the $q=0$ level, there is no more shares to sell,
having

$H(t,S,0)=0.$ (13)

  * •

If the price reaches a maximum (we have to fix a maximum price $S_{max}$ for
the numerical approach) the revenue is the amount of shares sold multiplied by
the fixed price. That is the inventory multiplied by the fixed maximum price.

$\displaystyle H(t,S_{max},q)=S_{max}q.$ (14)

  * •

We avoid the case of having any inventory $q>0$ at final time $T$, thus we
introduce a heavy penalty (minus infinity to be sure it never happens)

$H(T,S,q)=-\infty.$ (15)

By using central differences on $S$, forward differences on $q$ and backward
differences on $t$, we obtain the following discretized differential equation
that must be optimized

$\begin{split}\frac{H_{i,j}^{k+1}-H_{i,j}^{k}}{\Delta
t}+\frac{\sigma^{2}}{2}\left(\frac{H_{i+1,j}^{k}-2H_{i,j}^{k}+H_{i-1,j}^{k}}{(\Delta
s)^{2}}\right)-\nu\left(\frac{H_{i,j}^{k}-H_{i,j-1}^{k}}{\Delta q}\right)\\\
-g(\nu)\left(\frac{H_{i+1,j}^{k}-H_{i-1,j}^{k}}{2\Delta
s}\right)+\left(i\Delta s-\frac{1}{2}\Delta-f(\nu)\right)\nu=0,\end{split}$
(16)

which has approximation of orders $O(\Delta t),O(\Delta s^{2}),O(\Delta q)$.

By grouping similar terms and leaving the unknown terms (time $k\Delta t$ and
space $\geq j\Delta q$) on the same side of the equality we obtain

$AH_{i-1,j}^{k}+BH_{i,j}^{k}+CH_{i+1,j}^{k}=Z_{i,j}^{k*}$ (17)

where

$Z_{i,j}^{k*}=DH_{i,j}^{k+1}+EH_{i,j-1}^{k}+F_{i}(\nu)$
$F_{i}(\nu)=-\left(i\Delta s-\frac{1}{2}\Delta-f(\nu)\right)\nu$

and $A$,$B$,$C$,$D$,$E$ are the corresponding constants.

Since there is no diffusion in the $q$ variable the traditional $ADI$ scheme
for $2$ dimensions could lead to oscillations and instabilities. That is why
we prefer this approach, namely Directional $2$-dimensional finite
differences, where we use implicit solutions for vectors changing the variable
$S$ but having the variable $q$ constant, which additionally simplifies to
tridiagonal systems. More precisely, for each time $k\Delta t$ before expiry,
we begin by calculating a boundary condition at $q=0$, and once solved for
$q=(j-1)\Delta q$, we step to $q=j\Delta q$ and use the tridiagonal system in
(17) to solve for all the values of $S$, and iterate in ascending order over
$q$ up to $Q_{max}$. That is, on a uniform grid, we find the values of the
function $H$ for all values of $i$ while iterating in ascending order over $j$
and then stepping backwards through time in $k$, as shown in figure (1).

$(0,0)$$H_{1,1}^{k+1}$$H_{1,0}^{k}$$H_{2,0}^{k}$$H_{3,0}^{k}$$H_{4,0}^{k}$$H_{5,0}^{k}$$H_{0,1}^{k}$$H_{1,1}^{k}$$H_{2,1}^{k}$$H_{3,1}^{k}$$H_{4,1}^{k}$$H_{5,1}^{k}$$H_{6,1}^{k}$.$S=0$.$q=0$.$H(t,S_{max},q)$$t=T-k\Delta
t$$Q$$S$ Figure 1: This figure shows an example of the uniform grid at time
$t=T-k\Delta t$, where the brown highlighted line corresponds to the boundary
at $S=0$, the yellow to the boundary at $S=S_{max}$ and the green to the
boundary at $q=0$. The blue dots correspond to known values and the red dots
correspond to unknown values that will be found simultaneously at each step
through $q$ using the tridiagonal system (17)

####  3.2.2 Optimization

For the optimization process, since $f(\nu)$ is the TPI and $g(\nu)$ is the
PPI, we apply optimal policy iteration using first order condition (FOC) on
the control $\nu$, that is

$-g^{\prime}(\nu)\partial_{s}H-\partial_{q}H+S-\frac{1}{2}\Delta-\left(f^{\prime}(\nu)\nu+f(\nu)\right)=0$
(18)

Specifically:

  * •

For linear TPI and quadratic PPI, that is $f\left(\nu\right)=a_{1}\nu+a_{2}\
$, and $g(\nu)=c_{1}\nu^{2}+c_{2}\nu+c_{3}$ the FOC is

$-\left(2c_{1}\nu+c_{2}\right)\partial_{s}H-\partial_{q}H+S-\frac{1}{2}\Delta-\left(2a_{1}\nu+a_{2}\right)=0.$
(19)

  * •

For power TPI and power PPI, that is
$f\left(\nu\right)=r_{1}\nu^{r_{2}}+r_{3}\ $, and
$g(\nu)=p_{1}\nu^{p_{2}}+p_{3}$ the FOC is

$-p_{1}p_{2}\partial_{s}H\nu^{p_{2}-1}-\partial_{q}H+S-\frac{1}{2}\Delta-\left(r_{1}(r_{2}+1)\nu^{r_{2}}+r_{3}\right)=0.$
(20)

Note that linear cases are included in the power form just by taking $r_{2}=1$
in $f(\nu)$ and/or $p_{2}=1$ in $g(\nu)$.

##  4 Model calibration

This section presents the implementation of a method, based on Cartea and
Jaimungal, 2016a and Cartea et al.,  (2015), to estimate functional forms of
TPI and PPI. Additionally, to estimate the volatility, we use the realized
volatility, i.e,
$\hat{\sigma}=\sqrt{\sum_{t=1}^{T}\left(\ln(S_{t})-\ln(S_{t-1})\right)^{2}}$,
which is consistent with the model.

###  4.1 Temporary Price Impact

To find the functional form of the TPI, following the ideas from Cartea et
al.,  (2015), we use LOB data with $n$ ticks of depth taken at $N$ intervals
of $\tau$ seconds. In each of those $N$ intervals, we simulate the execution
of different sizes of Sell Market Orders, i.e. we liquidate the $m$ volumes
(As orders of different sizes) $Q=\\{Q_{1},Q_{2},...Q_{i},...Q_{m}\\}$ against
the existing LOB, walking the book when necessary. The LOB is represented, for
each available tick $j$, by a price $S_{j}$ and its corresponding number of
Limit Orders $Q_{S_{j}}$. These volumes are chosen in an increasing order up
to a significant maximum relative to the LOB’s total capacity (all posted
limit orders to buy). To find the execution price we define it as the weighted
average of the execution prices for each volume, in bid side, that is

$\displaystyle\hat{S}^{bid}_{i,t}=\frac{\sum_{j=1}^{N_{i}}S_{j}Q_{s_{j}}^{*}}{Q_{i}},\quad\text{
where }Q_{s_{j}}^{*}=\begin{cases}Q_{i}-\sum_{k=1}^{j-1}Q_{s_{k}}&\text{if
}Q_{i}\leq\sum_{k=1}^{j}Q_{s_{k}}\\\ Q_{s_{j}}&\text{else.}\end{cases},$

and $\displaystyle N_{i}=\min\left\\{j,\sum_{k=1}^{j}Q_{s_{k}}\geq
Q_{i}\right\\}$, and $0\leq t\leq N_{t}\tau$. Since we perform this procedure
every $\tau$ seconds, we consider not only volumes, but also liquidation rates
($\nu$) in our simulations, i.e. volumes per $\tau$ seconds. Thus, by letting
$S_{t}^{bid}$ be the best price on the bid side, the TPI at time $t$ due to
the liquidation rate $\displaystyle\nu_{i}=\frac{Q_{i}}{\tau}$ is
$\displaystyle TPI^{bid}_{i,t}=\hat{S}^{bid}_{i,t}-S^{bid}_{t}$. Finally, we
use the average of $TPI^{bid}_{i,t}$ over time, denoted by $TPI^{bid}_{i}$,
versus $\nu_{i}$ to fit the models $TPI^{bid}(\nu)=f(\nu)=a_{1}\nu+a_{2}$ and
$TPI^{bid}(\nu)=f(\nu)=r_{1}\nu^{r_{2}}+r_{3}$, via least squares.

###  4.2 Permanent Price Impact

To estimate a functional form of the permanent price impact, we use the same
idea as for the temporary impact, i.e. going through the data every $\tau$
seconds and at each of those moments executing simulated Market Sell Orders of
different volume sizes and therefore different trading liquidation rates. Then
we capture the impact of each trading rate, in this case, on the midprice
$\displaystyle S^{mid}_{t}=\frac{S^{bid}_{t}+S^{ask}_{t}}{2}$, that is, the
difference between the mean price existing before the liquidation of $Q_{i}$
shares in $\tau$ seconds and the mean price remaining after the liquidation.
Using least squares, we estimate the parameters of the models
$PPI^{bid}(\nu)=g(\nu)=b_{1}\nu+b_{2}$ and
$PPI^{bid}(\nu)=g(\nu)=p_{1}\nu^{p_{2}}+p_{3}$.

##  5 Results

###  5.1 Model parameters calibration

We use LOB intraday data with $n=100$ ticks of depth and a frequency of
$\tau=5$ seconds for $N_{t}=1440$ (two hours) taken from Binance API for BNB
666BNB is the official cryptocurrency of the Binance exchange and powers the
BNB Chain ecosystem on February 6, 2022.

A first aspect to consider in the calibration is to capture, from the LOB, the
total available volume at each moment that could be liquidated, i.e., for each
moment add up the available volumes (number of shares per tick) of the 100
ticks. This information is very important in the calibration of the functional
form of both the TPI and the PPI, because if the amount of shares to be
liquidated over a time lapse $\tau$ is small compared to the available volumes
in the LOB, the impact will be underestimated, since it would barely walk the
book. In this case, the calibration generates a power functional form
$(x^{\alpha}\text{ with }\alpha>1)$, we call this scenario underestimation. In
the opposite case, if the amount of shares to be liquidated is higher than
available volumes, the impact will be overestimated, since transactions have
to walk the book. Moreover, after a certain volume the full depth of the LOB
will be reached and therefore the impact will not grow more, this leads to a
power functional form $(x^{\alpha}\text{ with }0<\alpha<1)$, we call this
scenario overestimation. Finally, the intermediate case is found thinking of
liquidation rates, we see the effect on $\nu=Q/\tau$ since liquidation volumes
are considered in a time lapse $\tau$. Thus, for underestimation scenario, we
consider $\nu_{max}=50$, for overestimation scenario $\nu_{max}=7000$, and as
an intermediate value between the two previous cases, we use $\nu_{max}=1200$.
This intermediate scenario leads to a calibration of the functional forms of
TPI and PPI very close to linear function $(x^{\alpha}\text{ with }\alpha=1)$,
we call this scenario average-estimation.

The underestimated scenario is shown in figure 2, where subfigure (2(a)) shows
the plot of trading rates $\nu_{i}$ vs $TPI^{bid}_{i}$, for volumes $1\leq
i\leq m$, as well as the plots of the linear and power calibrated models
$TPI^{bid}(\nu)=f(\nu)=a_{1}\nu+a_{2}$ and
$TPI^{bid}(\nu)=f(\nu)=r_{1}\nu^{r_{2}}+r_{3}$. Subfigure (2(b)) shows plots
of trading rates $\nu_{i}$ vs. $PPI^{bid}_{i}$ and the calibrated PPI models
$PPI^{bid}(\nu)=g(\nu)=b_{1}\nu+b_{2}$ and
$PPI^{bid}(\nu)=g(\nu)=p_{1}\nu^{p_{2}}+p_{3}$. We can see that the shape of
the average TPI and PPI function can be calibrated more closely through a
power function than through a linear function. Due to the above, we consider,
to find optimal policies in this scenario, power $(x^{\alpha}\text{ with
}\alpha>0\text{ and }\alpha\neq 1)$ and linear $(x^{\alpha}\text{ with
}\alpha=1)$ impact forms.

(a) TPI

(b) PPI

Figure 2: Underestimation of power and linear functional forms of TPI and PPI
with $\nu_{max}=50$, $m=50$

In figure 3, we have the overestimated scenario and specifically note the
changes in the concavity of the calibrations.

(a) TPI

(b) PPI

Figure 3: Overestimation of power and linear functional forms of TPI and PPI
with $\nu_{max}=7000$,$m=50$

Additionally, the average-estimation scenario is in figure 4.

(a) TPI

(b) PPI

Figure 4: Average-estimation of power and linear functional forms of TPI and
PPI with $\nu_{max}=1200$,$m=50$

To measure which of the two functional forms, linear or power, fits best the
average TPI or PPI on the different scenarios, we use the squared residuals,
defined as $\hat{\epsilon}^{2}_{i,t}=\left(y_{i}-\hat{y}_{i}\right)^{2}_{t}$
and the corresponding R-squared, where $y_{i}$ is $TPI^{bid}_{i}$ or
$PPI^{bid}_{i}$ and $\hat{y_{i}}$ is $TPI^{bid}(\nu_{i})$ or
$PPI^{bid}(\nu_{i})$ (respectively). In appendix C, we show these results and
we observe that, in all cases, both the total sum of the residuals and the
average of the residuals are lower for the power functional form than for the
linear form. Likewise, the $R^{2}$ is higher in all cases for the power
functional form than for the linear form. This implies a better fit in all
cases for the power functional form, even in the average scenario where the
differences are minimal. We also obtain an average spread of $\Delta=0.100069$
and an estimated constant volatility of $\sigma=0.009388$.

###  5.2 Optimal policies

This section refers to the algorithm optimizing the trading rate, that is the
numerics for solving the $\mathcal{H}\mathcal{J}\mathcal{B}$ PDE and using the
optimal policies in (18) and (20).

####  5.2.1 Error and Convergence analysis

To verify our numerical algorithm, we perform some tests, first we calculate
the absolute error (when feasible) comparing the analytical and numerical
value function $H$ at $t=0$, that is for linear and quadratic PPI, obtaining
an accuracy of $10^{-5}$ and $10^{-4}$ respectively.

Additionally, we were able to confirm the convergence rate of the numerical
method for each of the independent variables of the value function $H$, i.e.
inventory $q$, price $s$ and time $t$, finding that the convergence rate is
$1$ for $q$ as expected. As for the price $s$, convergence rate is very close
to the theoretical value $2$ for all points of the grid, likewise, the
convergence rate for time $t$ is around its theoretical value $1$.

Finally, to verify that there are no significant differences between the
theoretical vs numerical optimal policies $\nu$, we applied the Chi-square
goodness of fit test between numerical and analytical policies for the linear
and quadratic PPI. We verify there are no significant differences between
these policies, i.e. in essence the policy found through the numerical method
is the same as the one found analytically. After performing these tests, we
are confident that our algorithm is correct and is working properly.

####  5.2.2 Analysis of optimal policies

Next, we show the optimal policies obtained using the numerical approach in
section 3.2, and we explore the results of using the calibrated linear and
power TPI and PPI for the scenarios under, over and average estimation.

Figure (6(a)) shows the plot of optimal policies along with a path for the
policy and inventory beginning at $0.5$ with power TPI and PPI in
underestimation scenario. To construct the path following the optimal policy
and the optimal inventory, starting at $q_{0}$ inventory, with the
corresponding optimal trading rate $\nu(0,q_{0})$, we find
$q_{1}=q_{0}-\nu(0,q_{0})\Delta t$. If $q_{1}$ lies on the grid, the procedure
continues, otherwise the respective value of $\nu(1,q_{1})$ is interpolated.
We repeat the procedure using $q_{k}=q_{k-1}-\nu(k-1,q_{k-1})\Delta t$ from
$t=0$ to $t=T=1$. Figure (6(b)) shows the corresponding 2D graph of the policy
path. Note that the figure (6(a)) is a plot of $\nu(t,q)$, i.e. not depending
on $S$, because we verify that optimal policies are independent of the price,
as we can see in figure (5).

(a) 3D optimal policy for fixed inventory $q=0.5$

(b) 3D optimal policy for fixed time $t=0.5$

Figure 5: Optimal policies for fixed inventory and time respectively. We
observe that the policy is price-independent in all scenarios and under the
different combinations of the TPI and PPI functions. The graphs show two
examples of this situation.

By changing the starting point of the inventory to $0.2$, figures (7(a)) and
(7(b)) show 3D optimal policy, optimal policy path, optimal inventory path and
corresponding 2D optimal policy path respectively. We observe some changes in
the liquidation policy, since starting with inventory $0.5$ implies
liquidation at rates close to $0$ after approximately time $t=0.9$, while when
starting with inventory $0.2$ this low rates are reached around $t=0.6$. One
possible interpretation of this situation is that, in an underestimation
scenario, a small impact causes the agent to liquidate early without much
impact and avoid the risk of exposure. We opted to show only these graphs
because we think are interesting since power TPI and PPI have best fit.
However, as the functional forms of TPI and PPI change, also does the policy.

(a) 3D optimal policy, optimal policy path and optimal inventory path

(b) 2D optimal policy path

Figure 6: Optimal policy for constant $S$ along with a path for the policy and
inventory when it begin at $0.5$ with power TPI and PPI in underestimation
scenario

(a) 3D optimal policy, optimal policy path and optimal inventory path

(b) 2D optimal policy path

Figure 7: Optimal policy for constant $S$ along with a path for the policy and
inventory when it begin at $0.2$ with power TPI and PPI in underestimation
scenario

Using parameters $T=1$, $S_{max}=300$, $Q_{max}=1$,$N_{t}=360$, $N_{s}=10$,
$N_{q}=100$, $\Upsilon=0.5$, the calibrated coefficients showed in appendix D
and using the mid point price $S=150$ as a constant, we plot the optimal
policy paths. In figure 8, we observe that if PPI is linear, the policy path
suggests to liquidate more slowly near initial time, even zero, to increase
the liquidation rate to a maximum near time $T$, in which the inventory is
zero, then the trading rate drop to zero. However, when PPI is power, the
policy starts liquidating quickly and decreasing the trading rate over time,
except towards the end of the liquidation period, because having positive
inventory causes the trading rate to increase slightly until all inventory is
liquidated. This may be caused by the fear of the heavy penalty imposed in the
boundary to avoid positive inventory at final time $T$.

We notice a huge difference between the policies, for different TPI and PPI
functions, and notably when $g$ changes from linear to power, because the
policy changes from increasing to decreasing. This may be caused by the
concavity induced in the underestimation scenario, since not walking the book
implies a faster execution and as there are less shares to liquidate the
trading rate decreases. However, this only happens when PPI is power, the
opposite happens when PPI is linear, which might suggest that in this scenario
linear PPI has a greater impact on the optimal policy, than power PPI.

Next we find the optimal policy paths for over and average-estimation of
linear and power TPI and PPI. Figures 9 and 10 show a comparison of the
optimal policy paths for different TPI and PPI. Interestingly, in the average-
estimation scenario we observe, that the optimal policy paths are similar if
PPI is power, regardless of the functional form of TPI, begins by not
liquidating and then increases until a certain point at which it drops to
zero.

For policies where PPI is linear, we observe that they start with a high
liquidation rate, decreasing over time and quickly reach a zero trading rate.
This is the opposite of the situation observed in the underestimation
scenario. This suggests that in the average-estimation scenario, power PPI has
a greater impact than linear PPI.

Similar situation is presented in the overestimation scenario, in figure 10,
only that in such scenario we observe that policy paths that liquidate early
(Linear PPI) do so more quickly than those liquidating late (power PPI).
Moreover, by having even a small inventory, the penalty imposed towards the
end of the time period causes a very large liquidation rate, which happens
when PPI is power.

According to the above, the optimal policy obtained is highly sensitive to the
model calibration process and therefore the optimal policy depends on the
calibration scenario, the functional form of both TPI and PPI and the LOB
data, since model coefficients depend entirely on this information.

Figure 8: Graphical comparison of optimal policies from underestimated TPI and
PPI, with $\Upsilon=0.5$ Figure 9: Graphical comparison of optimal policies
from average-estimated TPI and PPI, with $\Upsilon=0.5$ Figure 10: Graphical
comparison of optimal policies from overestimated TPI and PPI, with
$\Upsilon=0.5$

###  5.3 Simulation

This section presents simulations of the liquidation of $\Upsilon$ shares of
BNB using the optimal policies calculated from the previous section. The
simulation consists of taking an initial inventory $\Upsilon$ and liquidating
that inventory over time $[0,T]$ walking the book when needed, following each
of the policy paths and then compare the behavior of the numerical optimal
strategies versus the naive strategy (NS), which consists of liquidating all
cryptocurrencies in a single prespecified time.

To refer to the numerical optimal strategies we use the following notation:
the first letter refers to the estimation scenario, i.e. under (U), over(O)
and average(A), the second letter refers to TPI and therefore will always be
T, the third letter refers to the functional form of TPI, i.e. linear (L) or
power (P), the fourth letter will always be P in reference to PPI and the
fifth and last letter refers to the functional type of PPI, i.e. L or P. For
example, the optimal policy obtained in the underestimation scenario with
linear TPI and power PPI is represented by the letters UTLPP. With these
references, we use the policy paths starting with inventory $0.5$ for UTLPL,
UTLPP, UTPPL, UTPPP, OTLPL, ATLPL, ATLPP, ATPPL, ATPPP, OTLPP, OTPPL and
OTPPP.

In these simulations we assume, in the sense of Cartea and Jaimungal, 2016b ,
that the order book is resilient, i.e. when executing a MO, at the next moment
in time the book has again cleared that order and therefore the order book is
taken as it has been downloaded, without any modification. We begin showing in
table (1) the cumulative revenue obtained by simulating the execution of each
strategy as a percentage of NS, when we liquidate $\Upsilon=4000$ BNB
cryptocurrencies 777we choose a value of $4000$ to ensure that this inventory
is liquidatable at a single point in time in $1800$ seconds at intervals of
$5$ seconds.

Strategy | Accumulated
---|---
| Revenue (%)
UTLPL | 1.010922
UTLPP | 1.009505
UTPPL | 1.011206
UTPPP | 1.008875
ATLPL | 1.008409
ATLPP | 1.011454
ATPPL | 1.008456
ATPPP | 1.011473
OTLPL | 1.006268
OTLPP | 1.010698
OTPPL | 1.007157
OTPPP | 1.011098
NS | 1.000000
Table 1: Cumulative revenue as a percentage of NS

From the results in table 1, we observe that interestingly the policy paths
that start with a slow trading rate are the ones that earn the highest
revenue. The policy path that generated the highest revenue was ATPPP with a
cumulative revenue of $1.1473\%$ above NS and in second place was ATLPP with a
revenue of $1.1454\%$ above NS. This, although surprising, may be justified by
the fact that the number of shares to liquidate gives an average estimation
when compared to the LOB. In average-estimation scenario, these policy paths
are the ones that start liquidating slowly, even zero, increasing the trading
rate until a maximum point and when the inventory is zero, dropping to zero.
Similarly, in the overestimation and underestimation scenarios, the policy
paths that obtained the highest cumulative revenue were those that started
with a low trading rate.

Likewise, regardless of the scenario, the policy paths that obtained the
lowest cumulative revenue were those that liquidated faster at the beginning,
such is the case of OTLPL which was the one that obtained the lowest revenue
with a $0.6268\%$ above NS and OTPPL, which cumulative revenue was $0.7157\%$
above NS.

This situation can be described in terms of optimal inventory. If an optimal
policy path starts by liquidating slowly and then increase the liquidation
rate, the optimal inventory path will decrease slowly at first and then
decrease more rapidly, i.e., it is a concave function.

In the opposite case, when the policy path liquidates quickly at the beginning
and then decreases the liquidation rate, the inventory path will decrease
quickly at the beginning and then decrease more slowly, i.e., a convex
function.

Then, the optimal inventory path is a non-increasing function, which when
concave describes an optimal policy path that, as shown before, generates
higher cumulative revenue than convex optimal inventory functions, see figure
(11).

Figure 11: Optimal inventory paths for which, when concave, such as ATPPP and
ATLPP, generate higher cumulative revenues than convex paths, such as OTPPL
and OTLPL.

As a way of validating this behavior of the strategies we vary the inventory
$\Upsilon$, that is we simulate the liquidation of different inventories
starting at $500$ and ending at $4000$ with increments of $500$. The main
result of this simulation is shown in figure (12), in which we observe that as
the size of the inventory to be liquidated increases, the percentage of
cumulative revenue over NS obtained by optimal strategies also increases. We
also note that the strategies with highest revenue are always the same, i.e.
ATPPP and ATLPP and the under-performing strategies are also the same, i.e.
OTPPL and OTLPL. This is an indication that the performance of the strategy is
independent of $\Upsilon$.

Figure 12: Comparison of the performance of the optimal strategies through
accumulated revenue as a percentage of NS when inventory varies from $500$ to
$4000$ with steps of $500$.

To confirm or reject the empirical hypothesis that optimal policy paths that
start with a slow liquidation rate and increase as time passes are the ones
generating higher cumulative revenue and that, on the contrary, policy paths
that start with a high trading rate and decrease as time passes are the ones
generating less cumulative revenue, we consider that the optimal inventory
paths can be modeled by the non-increasing monotonic function

$\begin{array}[]{lccc}q:&\left\\{t\in\mathbb{Z}:0\leq t\leq
T\right\\}&\rightarrow&[0,\Upsilon]\\\
&t&\rightarrow&d_{1}t^{d_{2}}+d_{3}\end{array}$ (21)

subject to $q(0)=\Upsilon$ and $q(T)=0$, then

$q(t)=\frac{-\Upsilon}{T^{d_{2}}}t^{d_{2}}+\Upsilon$ (22)

For the parameters $\Upsilon=4000$, $T=1800$ and taking the exponent $d_{2}$
in $\left[0.1,0.5,1,5,10\right]$ to obtain changes in the concavity of the
function $q(t)$, we build simulations of optimal inventory paths, see figure
(13(a)). With each of these inventory paths, we simulate the execution of the
corresponding optimal policy path, thus obtaining its cumulative revenue as a
percentage of NS, characterized by $d_{2}$, see figure (13(b)).

(a) Simulated optimal inventory paths from equation
$q(t)=\frac{-4000}{1800^{d_{2}}}t^{d_{2}}+4000$

(b) Cumulative revenue from simulated optimal inventory paths

Figure 13: Simulation of optimal inventory paths and corresponding cumulative
revenue as a percentage of NS.

We observe that when the exponent $d_{2}$ is less than 1, i.e, $q(t)$ is
convex, the cumulative revenue as a percentage of NS obtained is less than the
other strategies. It seems that by increasing the value of the exponent
$d_{2}$, the cumulative revenue also increases. However, if the value of the
exponent $d_{2}$ is very close to zero or very large, these inventory paths
begin to approximate those of the NS strategy and therefore the corresponding
cumulative revenue obtained decrease. Figure (14) shows the cumulative revenue
as a percentage of NS obtained by simulating optimal inventories paths when
the exponent $d_{2}$ varies from $0.01$ to $100$ with increments of $0.1$. We
can observe that for values of $d_{2}$ close to zero, the cumulative revenue
is close to that obtained by NS and that for large values of the exponent
$d_{2}$ the cumulative revenue decreases. This situation leads us to conclude
that in an interval around $1$ for the exponent $d_{2}$, the concave optimal
inventories paths generate a higher cumulative revenue than convex optimal
inventories paths.

Figure 14: Cumulative revenue from simulated optimal inventories paths for
$d_{2}$ from $0.01$ to $100$ with steps of $0.1$, liquidating $\Upsilon=4000$.

In terms of optimal policies, those that start with a slow trading rate and
increase over time, as opposed to those that start liquidating more quickly
and decrease over time, generate higher cumulative revenue. Likewise, we
emphasize that this behavior is only observed if such policy paths do not
reach the extreme of resembling the NS, since in that case their cumulative
revenue decreases.

##  6 Discussion and conclusions

In this work we study the problem of optimal liquidation by analyzing the
consequences of having different functional forms for temporary impact,
specifically linear and power, and for permanent impact, specifically linear,
quadratic and more general power. We apply this theory in the cryptocurrency
market given its importance due to its evident growth, high volatility and its
articulation in the world of Decentralized Finance (DeFi) and secondly due to
the ease of accessing historical and LOB information. Therefore, we use LOB
data for BNB cryptocurrency, information downloaded from the Binance API.

When using LOB data, the depth and therefore the volumes available at the LOB
play an important role in the estimation and calibration process of both the
functional forms of TPI and PPI. That is, when simulating to liquidate
different inventories, if the maximum inventory to be liquidated is relatively
small to the inventory available in the data, the impact will be small at the
beginning of the simulation, but will increase as the size of the inventories
increases, which configures a power functional form ($x^{a}$ with
$a>1$).Additionally, if the inventory is larger than the available inventory
at the LOB, then at the beginning of the simulation there is a greater impact
than towards the end of the simulation, because the available inventory will
quickly be reached and therefore the impact will decrease, which configures a
power functional form ($x^{a}$, but with $0<a<1$). Likewise, for certain
values of inventory in the intermediate of these two extreme cases, the
functional form of the impacts tends to linearize ($x^{a}$ with $a=1$). Future
research aims to find specific characteristics in which each of the scenarios
is configured, as well as possible analytical models that allow identifying,
for example, points at which the impacts are linearized and differences
between these.

We observe that the shape of the optimal policy paths can change depending on
the starting inventory, the functional forms of TPI and PPI, and the values of
the calibrated parameters. However, the functional description of the policy
paths does not seem obvious, so we carry out this characterization through the
concavity of the optimal inventory paths. Leading us to acknowledge the
importance of the ratio between the initial amount of shares to be liquidated
and the availability of shares in the LOB when characterizing the performance
of optimal policies, and that these do not depend strongly on the functional
form of TPI and PPI.

Since this is a problem in which an amount $\Upsilon$ of cryptocurrencies is
to be liquidated in a pre-established period of time $T$, the optimal
inventory path can be described through a non-increasing monotonic function
$\displaystyle q(t)=\frac{-\Upsilon}{T^{d_{2}}}t^{d_{2}}+\Upsilon$, which
describes a large number of optimal policy paths. Under these conditions, we
observe that concave optimal inventory paths, i.e., $d_{2}>1$, are the ones
that generate the greatest cumulative revenue, under the condition that the
concavity is not so dramatic as the NS strategy.

This leads us to characterize the optimal winning policy paths as those whose
trading rate is low at the beginning and increases over time. In contrast to
those policy paths whose trading rate is high at the beginning and decreases
over time, which can be considered as losers, are characterized by convex
optimal inventory paths, that is, $0<d_{2}<1$.

Further research tends to examine this phenomenon in greater detail,
identifying, for example, the exact conditions and points at which the
described behavior changes, or, based on families of optimal policies,
identifying the “most” optimal among them, the functional forms of these
families of optimal policies, among other options.

The scope of this study is limited, initially by the LOB data available, by
the liquidity of the asset to be liquidated since high liquidity is required
so that the inventories can be liquidated, in this case in time horizons of
seconds. Other limitations of this study are the model used does not take into
consideration other variables that could affect the optimal liquidation policy
such as volatility variability, liquidating using limit orders, new market and
limit order flows by other agents, the liquidation dynamics that may imply
jumps due to ticks, among others. Therefore, considering several additional
variables to the problem, different models, different characteristics of
financial assets and markets in which they are traded are part of the future
research agenda.

## References

  * Alfonsi and Schied, (2010) Alfonsi, A. F. A. and Schied, A. (2010).  Optimal execution strategies in limit order books with general shape functions.  Quantitative Finance, Taylor.
  * Almgren and Chriss, (2000) Almgren, R. and Chriss, N. (2000).  Optimal execution of portfolio transactions.  J. Risk, 3:5–39.
  * Almgren et al., (2005) Almgren, R., Thum, C., Hauptmann, E., and Li, H. (2005).  Direct estimation of equity market impact.
  * Barger and Lorig, (2018) Barger, W. and Lorig, M. (2018).  Optimal liquidation under stochastic price impact.  International Journal of Theoretical and Applied Finance.
  * Barles and Souganidis, (1991) Barles, G. and Souganidis, P. E. (1991).  Convergence of approximation schemes for fully nonlinear second order equations.  Asymptotic Analysis, 4(3):271–283.
  * Bellman and Dreyfus, (1962) Bellman, R. and Dreyfus, S. (1962).  Applied dynamic programming.  A report kprepared for United States Air Force project RAND.
  * Bershova and Rakhlin, (2013) Bershova, N. and Rakhlin, D. (2013).  The non-linear market impact of large trades: Evidence from buy-side order flow.  Quantitative Finance, 13(11):1759–1778.
  * Blanco Encinosa, (2021) Blanco Encinosa, L. J. (2021).  Criptomonedas. Breve análisis desde la perspectiva económica y financiera.  Cofin Habana, 15.
  * Brunovský et al., (2018) Brunovský, P., Černý, A., and Komadel, J. (2018).  Optimal trade execution under endogenous pressure to liquidate: Theory and numerical solutions.  European Journal of Operational Research, 264(3):1159–1171.
  * (10) Cartea, A. and Jaimungal, S. (2016a).  A closed-form execution strategy to target.  SIAM Journal on Financial Mathematics, 7(1):760–785.
  * (11) Cartea, A. and Jaimungal, S. (2016b).  Incorporating order-flow into optimal execution.  Mathematics and Financial Economics, 10(3):339–364.
  * Cartea et al., (2015) Cartea, A., Jaimungal, S., and Penalva, J. (2015).  Algorithmic and high-frequency trading.  Cambridge University Press.
  * Curato et al., (2017) Curato, G., Gatheral, J., and Lillo, F. (2017).  Optimal execution with nonlinear transient market impact.  Quantitative Finance, 17(1):41–54.
  * Duffy, (2006) Duffy, D. J. (2006).  Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach.  John Wiley.
  * Gatheral, (2010) Gatheral, J. (2010).  No-dynamic-arbitrage and market impact.  Quantitative Finance, 10(7):749–759.
  * Gueant, (2014) Gueant, O. (2014).  Permanent market impact can be nonlinear.  Preprint.
  * Huberman and Stanzl, (2004) Huberman, G. and Stanzl, W. (2004).  Price manipulation and quasi-arbitrage.  Econometrica, 72(4):1247–1275.
  * Lillo et al., (2003) Lillo, F., Farmer, J. D., and Mantegna, R. N. (2003).  Master curve for price-impact function.  Nature, 421(1):129–130.
  * Peaceman and Ratchford Jr, (1955) Peaceman, D. and Ratchford Jr, H. H. (1955).  The numerical solution of parabolic and elliptic differential equations.  Journal of the Society for Industrial and Applied Mathematics, 3(1):28–41.
  * Pham, (2009) Pham, H. (2009).  Continuous-time Stochastic Control and Optimization with Financial Applications.  Springer.
  * Piñeiro-Chousa et al., (2022) Piñeiro-Chousa, J., Ángeles López-Cabarcos, M., Sevic, A., and González-López, I. (2022).  A preliminary assessment of the performance of defi cryptocurrencies in relation to other financial assets, volatility, and user-generated content.  Technological Forecasting and Social Change, 181:121740.
  * Subramanian, (2008) Subramanian, A. (2008).  Optimal liquidation by a large investor.  SIAM Journal of Applied Mathematics, 68:1168–1201.
  * Tavella and Randall, (2000) Tavella, D. and Randall, C. (2000).  Pricing Financial Instruments: The Finite Difference Method.  John Wiley.
  * Thomas, (1995) Thomas, J. (1995).  Numerical Partial Differential Equations: Finite Difference Methods.  Springer New York, NY.
  * Toth et al., (2016) Toth, B., Eisler, Z., and Bouchaud, J. P. (2016).  The square-root impact law also holds for option markets.  Wilmott, (85):70–73.

##  Appendix A Proof of Lemma 1

In this section we assume that $f(\nu)=a_{1}\nu+a_{2}$ and $g(\nu)=b\nu$. The
objective is to solve the $\mathcal{H}\mathcal{J}\mathcal{B}$ differential
equation (5) to find the policy that allows the optimal liquidation of $q_{t}$
shares in a period of time $[t,T)$. Thus we have to solve

${\partial}_{t}H+\frac{1}{2}{\sigma}^{2}{\partial}_{ss}H+\sup_{\nu\in\mathcal{A}}\left\\{-b\nu\partial_{s}H-\nu\partial_{q}H+\left(S-\frac{1}{2}\Delta-
a_{1}\nu-a_{2}\right)\nu\right\\}=0$ (23)

Using the first order conditions on $\nu$ we obtain:

$\partial_{t}H+\frac{1}{2}{\sigma}^{2}{\partial}_{ss}H+\frac{\left(S-\frac{1}{2}\Delta-
a_{2}-\partial_{q}H-b\partial_{s}H\right)^{2}}{4a_{1}}=0$ (24)

Using the educated guess

$H\left(t,S,q\right)=q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\frac{b}{2}q^{2}+h(t,q)$ (25)

and replacing in the differential equation (24) we simplify to

${\partial}_{t}h+\frac{1}{4a_{1}}{\left({\partial}_{q}h\right)}^{2}=0$ (26)

thus finding the following solutions

$H(t,s,q)=q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\left(\frac{b}{2}+\frac{a_{1}}{(T-t)}\right)q^{2}$ (27)
$\nu^{*}_{t}=\frac{q_{t}}{T-t}$ (28)

##  Appendix B Proof of Lemma 2

In this section we assume that $f(\nu)=a_{1}\nu+a_{2}$ and
$g(\nu)=c_{1}\nu^{2}+c_{2}\nu$. The objective is to solve the
$\mathcal{H}\mathcal{J}\mathcal{B}$ differential equation (5) to find the
optimal policy for liquidating $q_{t}$ shares in a time period $[t,T)$. Thus
we have to solve

${\partial}_{t}H+\frac{1}{2}{\sigma}^{2}{\partial}_{ss}H+\sup_{\nu\in\mathcal{A}}\left\\{\left(-c_{1}\nu^{2}-c_{2}\nu\right)\partial_{s}H-\nu\partial_{q}H+\left(S-\frac{1}{2}\Delta-
a_{1}\nu-a_{2}\right)\nu\right\\}=0$ (29)

Using the first order conditions on $\nu$ we obtain:

$\partial_{t}H+\frac{1}{2}{\sigma}^{2}{\partial}_{ss}H+\frac{\left(S-\frac{1}{2}\Delta-
a_{2}-\partial_{q}H-c_{2}\partial_{s}H\right)^{2}}{4\left(c_{1}{\partial}_{s}H+a_{1}\right)}=0$
(30)

Using the educated guess

$H\left(t,S,q\right)=q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\frac{c_{2}}{2}q^{2}+h(t,q)$ (31)

and replacing in the differential equation (30) we simplify to

${\partial}_{t}h+\frac{\left({\partial}_{q}h\right)^{2}}{4\left(c_{1}q+a_{1}\right)}=0$
(32)

thus finding the following solutions

$H\left(t,S,q\right)=\left\\{\begin{array}[]{lcc}q\left(S-\frac{1}{2}\Delta-
a_{2}\right)-\frac{c_{2}}{2}q^{2}-\frac{4\left(c_{1}q+a_{1}\right)^{3}}{9c_{1}^{2}(T-t)}&if&q\neq
0\\\ \\\ 0&if&q=0\end{array}\right.$ (33)
$\nu_{t}^{*}=\frac{2\left(c_{1}q_{t}+a_{1}\right)}{3c_{1}\left(T-t\right)}$
(34)

##  Appendix C Comparison squared residuals and $R^{2}$

Comparison squared residuals $(\hat{\epsilon}^{2})$ |  R-squared $(R^{2})$
---|---
impact | model | total sum | mean | std dev |
Underestimate TPI (UTPI) | Linear | 0.0000590 | 0.0000012 | 0.0000012 | 0.990627
Underestimate TPI (UTPI) | Power | 0.0000162 | 0.0000003 | 0.0000006 | 0.997424
Underestimate PPI (UPPI) | Linear | 0.0001789 | 0.0000037 | 0.0000029 | 0.980279
Underestimate PPI (UPPI) | Power | 0.0000171 | 0.0000003 | 0.0000004 | 0.998111
Overestimate TPI (OTPI) | Linear | 2.4095625 | 0.0491747 | 0.0738702 | 0.980702
Overestimate TPI (OTPI) | Power | 1.4869937 | 0.0303468 | 0.0374643 | 0.988091
Overestimate PPI (OPPI) | Linear | 5.6188366 | 0.1146701 | 0.1361679 | 0.955715
Overestimate PPI (OPPI) | Power | 3.8335537 | 0.0782358 | 0.0732665 | 0.969786
Average-estimate TPI (ATPI) | Linear | 0.0010416 | 0.0000213 | 0.0000237 | 0.999803
Average-estimate TPI (ATPI) | Power | 0.0010396 | 0.0000212 | 0.0000230 | 0.999804
Average-estimate PPI (APPI) | Linear | 0.0096960 | 0.0001979 | 0.0001787 | 0.998061
Average-estimate PPI (APPI) | Power | 0.0056630 | 0.0001156 | 0.0001251 | 0.998867
Table 2: Comparison of linear and power model fit through squared residuals
and $R^{2}$ for under, over and average estimated TPI and PPI

##  Appendix D Calibrated coeficients

Since the differential equation (5) incorporated the direction of the
negotiations, i.e. liquidation, through the negative sign for both TPI and
PPI, these coefficients must be considered through their additive inverse,
except those representing exponents, i.e. $r_{2}$ and $p_{2}$.

Coefficients obtained for each calibrated model
---
impact | model | First parameter | Second parameter | Third parameter
UTPI | Linear | $a_{1}=-0.00079754$ | $a_{2}=-0.00066177$ | NA
UTPI | Power | $r_{1}=-2.44114118e-04$ | $r_{2}=1.29520174$ | $r_{3}=-3.77323856e-03$
UPPI | Linear | $b_{1}=-0.00095264$ | $b_{2}=0.00229332$ | NA
UPPI | Power | $p_{1}=-1.38745021e-04$ | $p_{2}=1.48472878$ | $p_{3}=-3.01507718e-03$
OTPI | Linear | $a_{1}=-0.00079843$ | $a_{2}=-0.22235319$ | NA
OTPI | Power | $r_{1}=-0.00538481$ | $r_{2}=0.78904313$ | $r_{3}=0.23917224$
OPPI | Linear | $b_{1}=-0.00079455$ | $b_{2}=-0.2597109$ | NA
OPPI | Power | $p_{1}=-0.00947337$ | $p_{2}=0.72704129$ | $p_{3}=0.38622943$
ATPI | Linear | $a_{1}=-0.00095984$ | $a_{2}=0.00209078$ | NA
ATPI | Power | $r_{1}=-0.0011318$ | $r_{2}=0.97757467$ | $r_{3}=0.01148375$
APPI | Linear | $b_{1}=-0.0009179$ | $b_{2}=-0.01720435$ | NA
APPI | Power | $p_{1}=-0.00200298$ | $p_{2}=0.89422254$ | $p_{3}=0.02904742$
Table 3: Calibrated oefficients