A Parameter Estimation Method Using Linear Response Statistics

Abstract

This paper presents a new parameter estimation method for Itô diffusions such that the resulting model predicts the equilibrium statistics as well as the sensitivities of the underlying system to external disturbances. Our formulation does not require the knowledge of the underlying system, however, we assume that the linear response statistics can be computed via the fluctuation–dissipation theory. The main idea is to fit the model to a finite set of “essential” statistics that is sufficient to approximate the linear response operators. In a series of test problems, we will show the consistency of the proposed method in the sense that if we apply it to estimate the parameters of the underlying model, then we must obtain the true parameters.

Appendices

Appendix 1: The Computation of $M_2$ for the Simplified Turbulence Model

Here, we provide the computational detail for obtaining (27). In particular,

$$\begin{aligned} \hat{k}''_{x}(0)=\frac{1}{\sigma _{eq}^{2}}\int _{\mathbb {R}^{3}} \tilde{\mathcal {L}}(\tilde{B}(x,x)+\tilde{L}x-\tilde{\Lambda }x)x^{\top } \tilde{p}_{eq}(x)\text {d}x. \end{aligned}$$

where $\tilde{B}(x,x)=(\tilde{B}_1x_2x_3,\tilde{B}_2x_1x_3,\tilde{B}_3 x_1x_2)$ and $\tilde{L}$, $\tilde{\Lambda }$ are $3\times 3$ matrices. We introduce the following notations.

• 1. Set $\tilde{D}(x)=\tilde{B}(x,x)+\tilde{L}x-\tilde{\Lambda }x$ which is the deterministic part of the system, and $\tilde{D}_{i}$ denotes the ith component of $\tilde{D}(x)$.

• 2. Set $\tilde{C}=(\tilde{L}-\tilde{\Lambda })$, then we have $\tilde{D}(x)=\tilde{B}(x,x)+\tilde{C}x$, and $\tilde{C}_i$ denotes the ith row of the matrix, which is a row vector.

With these notations, we have the Jacobian matrix of $\tilde{D}(x)$

$$\begin{aligned} \nabla \tilde{D}(x)=\nabla \tilde{B}(x,x)+\nabla (\tilde{C}x)= \left( \begin{array}{lll} 0 &{}\quad \tilde{B}_1 x_3 &{}\quad \tilde{B}_1 x_2 \\ \tilde{B}_2 x_3 &{}\quad 0 &{}\quad \tilde{B}_2 x_1 \\ \tilde{B}_3 x_2 &{}\quad \tilde{B}_3 x_1 &{}\quad 0 \end{array} \right) +\tilde{C}, \end{aligned}$$

and the generator $\tilde{\mathcal {L}}$ acting on f becomes,

$$\begin{aligned} \tilde{\mathcal {L}}f =\sum _{i=1}^{3}\tilde{D}_{i}\frac{\partial f}{\partial x_{i}}+\frac{\tilde{\sigma }^{2}}{2}\sum _{i,j=1}^{3} (\tilde{\Lambda })_{i,j}\frac{\partial ^{2}f}{\partial x_{i}\partial x_{j}}. \end{aligned}$$

In our case, $f=\tilde{D}(x)$. Since $\tilde{p}_{eq}\sim \mathcal {N}(0,\tilde{\sigma }_{eq}^{2})$, with $\tilde{\sigma }_{eq}^{2}=\frac{\tilde{\sigma }^{2}}{2}$, only the odd power terms in $\tilde{\mathcal {L}}\tilde{D}(x)$ contribute to $M_2$. Notice that the second order partial derivatives of $\tilde{D}(x)$ are constant, thus,

$$\begin{aligned} \hat{k}''_{x}(0)=\frac{1}{\sigma _{eq}^{2}}\int _{\mathbb {R}^{3}} (\tilde{\mathcal {L}}\tilde{D}(x))x^{\top }\tilde{p}_{eq}\text {d}x= \frac{1}{\sigma _{eq}^{2}}\int _{\mathbb {R}^{3}}\left( \sum _{i=1}^{3} \tilde{D}_{i}\frac{\partial \tilde{D}}{\partial x_{i}}\right) x^{\top }\tilde{p}_{eq}\text {d}x. \end{aligned}$$

For example, if $i=1$,

$$\begin{aligned} \tilde{D}_{1}\frac{\partial \tilde{D}}{\partial x_1}&=(\tilde{B}_1x_2x_3+\tilde{C}_1 x)[(0,\tilde{B}_2x_3,\tilde{B}_3x_2)^{\top }+\tilde{C}_1^{\top }]\\&=(0,\tilde{B}_1\tilde{B}_2x_2x_3^2,\tilde{B}_1 \tilde{B}_3x_2^2x_3)^{\top }+\tilde{B}_1x_2x_3\tilde{C}_1^{\top } +\tilde{C}_1x(0,\tilde{B}_2x_3,\tilde{B}_3x_2)^{\top }+\tilde{C}_1x\tilde{C}_1^{\top }. \end{aligned}$$

Notice that both $\tilde{B}_1x_2x_3\tilde{C}_1^{\top }$ and $\tilde{C}_1x(0,\tilde{B}_2x_3,\tilde{B}_3x_2)^{\top }$ are even power terms which means that they do not affect the value of $M_2$. Then we have,

$$\begin{aligned} \int _{\mathbb {R}^{3}}\tilde{D}_{1}\frac{\partial \tilde{D}}{\partial x_1}x^{\top }\tilde{p}_{eq}\text {d}x&= \int _{\mathbb {R}^{3}} \left( \begin{array}{c} 0 \\ \tilde{B}_1\tilde{B}_2 x_2x_3^2 \\ \tilde{B}_1\tilde{B}_3 x_2^2x_3 \end{array} \right) x^{\top }\tilde{p}_{eq}+\tilde{C}_{1}^{\top } \tilde{C}_{1}xx^{\top }\tilde{p}_{eq}\text {d}x \\&= \int _{\mathbb {R}^{3}} \left( \begin{array}{c c c} 0 &{}\quad 0 &{}\quad 0 \\ \tilde{B}_1\tilde{B}_2 x_1x_2x_3^2 &{}\quad \tilde{B}_1\tilde{B}_2x_2^2x_3^2 &{}\quad \tilde{B}_1\tilde{B}_2 x_2x_3^3\\ \tilde{B}_1\tilde{B}_3 x_1x_2^2x_3 &{}\quad \tilde{B}_1\tilde{B}_3 x_2^3x_3 &{}\quad \tilde{B}_1\tilde{B}_3 x_2^2x_3^2 \end{array} \right) \tilde{p}_{eq} \text {d}x\\&\quad +\,\tilde{C}_{1}^{\top }\tilde{C}_{1}\int _{\mathbb {R}^{3}} xx^{\top }\tilde{p}_{eq}\text {d}x \\&= \tilde{\sigma }_{eq}^4 \left( \begin{array}{c c c} 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \tilde{B}_1\tilde{B}_2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \tilde{B}_1\tilde{B}_3 \end{array} \right) +\tilde{\sigma }_{eq}^2\tilde{C}_{1}^{\top }\tilde{C}_{1}. \end{aligned}$$

Similarly, for $i=2$ and $i=3$, we have

$$\begin{aligned} \int _{\mathbb {R}^{3}}\tilde{D}_{2}\frac{\partial \tilde{D}}{\partial x_2}x^{\top }\tilde{p}_{eq}\text {d}x = \tilde{\sigma }_{eq}^4 \left( \begin{array}{c c c} \tilde{B}_1\tilde{B}_2 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \tilde{B}_2\tilde{B}_3 \end{array} \right) +\tilde{\sigma }_{eq}^2\tilde{C}_{2}^{\top }\tilde{C}_{2}, \\ \int _{\mathbb {R}^{3}}\tilde{D}_{3}\frac{\partial \tilde{D}}{\partial x_3}x^{\top }\tilde{p}_{eq}\text {d}x = \tilde{\sigma }_{eq}^4 \left( \begin{array}{c c c} \tilde{B}_1\tilde{B}_3 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \tilde{B}_2\tilde{B}_3 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{array} \right) +\tilde{\sigma }_{eq}^2\tilde{C}_{3}^{\top }\tilde{C}_{3}, \end{aligned}$$

respectively. Thus,

$$\begin{aligned} \hat{k}''_{x}(0)=\frac{\tilde{\sigma }_{eq}^{2}}{\sigma _{eq}^{2}} \left[ \tilde{\sigma }_{eq}^{2} \left( \begin{array}{c c c} \tilde{B}_1\tilde{B}_2+\tilde{B}_1\tilde{B}_3 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad \tilde{B}_1\tilde{B}_2+\tilde{B}_2\tilde{B}_3 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \tilde{B}_1\tilde{B}_3+\tilde{B}_2\tilde{B}_3 \end{array} \right) +\tilde{C}^{2} \right] . \end{aligned}$$

Furthermore, since $\sum _{i} \tilde{B}_i=0$, we arrive at,

$$\begin{aligned} \hat{k}''_{x}(0)=\frac{\tilde{\sigma }_{eq}^{2}}{\sigma _{eq}^{2}} [-\tilde{\sigma }_{eq}^{2} \text {diag}(\tilde{B}_1^2, \tilde{B}_2^2, \tilde{B}_3^2 )+\tilde{C}^2], \end{aligned}$$

which is the stated result.

Appendix 2: The Computation of $M_4$ and $M_5$ for the Langevin Dynamics Model

Here we show how to obtain the formulas for $\hat{k}^{(4)}_A(0)$ and $\hat{k}^{(5)}_A(0)$ for the Langevin dynamics model. Recall that the generator is given by,

$$\begin{aligned} \tilde{\mathcal {L}} = v\frac{\partial }{\partial x} + (-U'(x)-\tilde{\gamma }v) \frac{\partial }{\partial v} + \tilde{\gamma }k_B\tilde{T} \frac{\partial ^2}{\partial v^2}. \end{aligned}$$

Furthermore, we have the response operator,

$$\begin{aligned} \hat{k}_{A}(t) = \frac{1}{k_B T}\int _{\mathbb {R}^2}\int _{\mathbb {R}^2} v u \tilde{p}(x,v,t|y,u,0) \text {d}x\text {d}v\text {d}y \text {d}u. \end{aligned}$$

Direct calculations yield,

$$\begin{aligned} \hat{k}_{A}'(t)&= \frac{1}{k_B T}\int _{\mathbb {R}^2}\int _{\mathbb {R}^2} (-U'(x)-\tilde{\gamma }v)u \tilde{p}(x,v,t|y,u,0) \text {d}x\text {d}v\text {d}y\text {d}u,\\ \hat{k}_{A}''(t)&= \frac{1}{k_B T}\int _{\mathbb {R}^2}\int _{\mathbb {R}^2} (-vU''(x)+\tilde{\gamma }U'(x)+\tilde{\gamma }^{2} v) u \tilde{p}(x,v,t|y,u,0) \text {d}x\text {d}v\text {d}y \text {d}u. \end{aligned}$$

From fitting $M_0$, we have $T=\tilde{T}.$ Now we applying the generator $\tilde{\mathcal {L}}$ again and we obtain

$$\begin{aligned} \hat{k}'''_A&=\frac{1}{k_B T}\int (-v^2 U'''(x)+2\tilde{\gamma }v U''(x)+U'(x)U''(x)-\tilde{\gamma }^{3}U'(x)\\&\quad -\,\tilde{\gamma }^3v)u \tilde{p}(x,v,t|y,u,0)\text {d}x \text {d}v\text {d}y\text {d}u, \end{aligned}$$

which leads to $k'''_A(0)=2\tilde{\gamma }\mathbb {E}_{eq}(U''(x))-\tilde{\gamma }^3$. And applying the generator again we have

$$\begin{aligned} \hat{k}^{(4)}_A=\frac{1}{k_B T}\int f(x,v)u \tilde{p}(x,v,t|y,u,0)\text {d}x \text {d}v\text {d}y\text {d}u, \end{aligned}$$

where $f(x,v)=-v^3U^{(4)}(x)+4\tilde{\gamma }v^2U'''(x)+v(U''(x))^{2}+3vU'(x)U''' (x)-2\tilde{\gamma }k_B\tilde{T}U'''(x)-2\tilde{\gamma }U'(x)U''(x)-3\tilde{\gamma }^{2}vU''(x) +\tilde{\gamma }^{3}U'(x)+\tilde{\gamma }^{4}v$. This formula leads to

$$\begin{aligned} \hat{k}^{(4)}_A(0)&=-\frac{1}{k_B\tilde{T}}\mathbb {E}_{eq}(v^4) \mathbb {E}_{eq}U^{(4)}(x)+\mathbb {E}_{eq}((U''(x))^2) +3\mathbb {E}_{eq}(U'(x)U'''(x))\\&\quad -\,3\tilde{\gamma }^2\mathbb {E}_{eq}(U''(x)) +\tilde{\gamma }^4. \end{aligned}$$

For $\hat{k}_A^{(5)}(0)$ we have

$$\begin{aligned} \hat{k}_A^{(5)}(0)&=\frac{7\tilde{\gamma }}{k_BT}\mathbb {E}_{eq}(v^4) \mathbb {E}_{eq}U^{(4)}(x)-8\tilde{\gamma }k_BT\mathbb {E}_{eq}(U^{(4)} (x))-3\tilde{\gamma }\mathbb {E}_{eq}((U''(x))^2) \\&\quad -\,13\tilde{\gamma }\mathbb {E}_{eq}(U'(x)U'''(x)) +4\tilde{\gamma }^3\mathbb {E}_{eq}(U''(x))-\tilde{\gamma }^5. \end{aligned}$$

Thus these formulas require the fourth moment of v, $\mathbb {E}_{eq}(U'U''')$, $\mathbb {E}_{eq}(U^{(4)})$, $\mathbb {E}_{eq}((U'')^2)$ and $\mathbb {E}_{eq}(U'')$, in order to compute $\hat{k}_A^{(4)}(0)$ and $\hat{k}_A^{(5)}(0)$.

Notice $v\sim \mathcal {N}(0,k_{B}T)$, thus $\mathbb {E}_{eq}(v^4)=3(k_B T)^2$. Therefore, we can rewrite the two formula into

$$\begin{aligned} \hat{k}^{(4)}_A(0)&=-3k_BT\mathbb {E}_{eq}U^{(4)}(x)+\mathbb {E}_{eq} ((U''(x))^2)+3\mathbb {E}_{eq}(U'(x)U'''(x))\\&\quad -\,3\tilde{\gamma }^2\mathbb {E}_{eq}(U''(x))+\tilde{\gamma }^4, \\ \hat{k}_A^{(5)}(0)&=13\tilde{\gamma }k_BT\mathbb {E}_{eq}(U^{(4)}(x)) -3\tilde{\gamma }\mathbb {E}_{eq}((U''(x))^2)-13\tilde{\gamma } \mathbb {E}_{eq}(U'(x)U'''(x))\\&\quad +\,4\tilde{\gamma }^3\mathbb {E}_{eq}(U''(x)) -\tilde{\gamma }^5. \end{aligned}$$

These are the formulas implemented in producing the results in Fig. 2.