def sde(self):
<EOL>if self.approx_order is not None:<EOL><INDENT>N = self.approx_order<EOL><DEDENT>else:<EOL><INDENT>N = <NUM_LIT:10><EOL><DEDENT>roots_rounding_decimals = <NUM_LIT:6><EOL>fn = np.math.factorial(N)<EOL>p_lengthscale = float( self.lengthscale )<EOL>p_variance = float(self.variance)<EOL>kappa = <NUM_LIT:1.0>/<NUM_LIT>/p_lengthscale**<NUM_LIT:2><EOL>Qc = np.array( ((p_variance*np.sqrt(np.pi/kappa)*fn*(<NUM_LIT:4>*kappa)**N,),) )<EOL>eps = <NUM_LIT><EOL>if (float(Qc) > <NUM_LIT:1.0>/eps) or (float(Qc) < eps):<EOL><INDENT>warnings.warn("""<STR_LIT>""".format(float(Qc)) )<EOL><DEDENT>pp1 = np.zeros((<NUM_LIT:2>*N+<NUM_LIT:1>,)) <EOL>for n in range(<NUM_LIT:0>, N+<NUM_LIT:1>): <EOL><INDENT>pp1[<NUM_LIT:2>*(N-n)] = fn*(<NUM_LIT>*kappa)**(N-n)/np.math.factorial(n)*(-<NUM_LIT:1>)**n<EOL><DEDENT>pp = sp.poly1d(pp1)<EOL>roots = sp.roots(pp)<EOL>neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < <NUM_LIT:0>]<EOL>aa = sp.poly1d(neg_real_part_roots, r=True).coeffs<EOL>F = np.diag(np.ones((N-<NUM_LIT:1>,)),<NUM_LIT:1>)<EOL>F[-<NUM_LIT:1>,:] = -aa[-<NUM_LIT:1>:<NUM_LIT:0>:-<NUM_LIT:1>]<EOL>L= np.zeros((N,<NUM_LIT:1>))<EOL>L[N-<NUM_LIT:1>,<NUM_LIT:0>] = <NUM_LIT:1><EOL>H = np.zeros((<NUM_LIT:1>,N))<EOL>H[<NUM_LIT:0>,<NUM_LIT:0>] = <NUM_LIT:1><EOL>Pinf = lyap(F, -np.dot(L,np.dot( Qc[<NUM_LIT:0>,<NUM_LIT:0>],L.T)))<EOL>Pinf = <NUM_LIT:0.5>*(Pinf + Pinf.T)<EOL>dF    = np.empty([F.shape[<NUM_LIT:0>],F.shape[<NUM_LIT:1>],<NUM_LIT:2>])<EOL>dQc   = np.empty([Qc.shape[<NUM_LIT:0>],Qc.shape[<NUM_LIT:1>],<NUM_LIT:2>])<EOL>dPinf = np.empty([Pinf.shape[<NUM_LIT:0>],Pinf.shape[<NUM_LIT:1>],<NUM_LIT:2>])<EOL>dFvariance = np.zeros(F.shape)<EOL>dFlengthscale = np.zeros(F.shape)<EOL>dFlengthscale[-<NUM_LIT:1>,:] = -aa[-<NUM_LIT:1>:<NUM_LIT:0>:-<NUM_LIT:1>]/p_lengthscale * np.arange(-N,<NUM_LIT:0>,<NUM_LIT:1>)<EOL>dQcvariance = Qc/p_variance<EOL>dQclengthscale = np.array(( (p_variance*np.sqrt(<NUM_LIT:2>*np.pi)*fn*<NUM_LIT:2>**N*p_lengthscale**(-<NUM_LIT:2>*N)*(<NUM_LIT:1>-<NUM_LIT:2>*N),),))<EOL>dPinf_variance = Pinf/p_variance<EOL>lp = Pinf.shape[<NUM_LIT:0>]<EOL>coeff = np.arange(<NUM_LIT:1>,lp+<NUM_LIT:1>).reshape(lp,<NUM_LIT:1>) + np.arange(<NUM_LIT:1>,lp+<NUM_LIT:1>).reshape(<NUM_LIT:1>,lp) - <NUM_LIT:2><EOL>coeff[np.mod(coeff,<NUM_LIT:2>) != <NUM_LIT:0>] = <NUM_LIT:0><EOL>dPinf_lengthscale = -<NUM_LIT:1>/p_lengthscale*Pinf*coeff<EOL>dF[:,:,<NUM_LIT:0>]    = dFvariance<EOL>dF[:,:,<NUM_LIT:1>]    = dFlengthscale<EOL>dQc[:,:,<NUM_LIT:0>]   = dQcvariance<EOL>dQc[:,:,<NUM_LIT:1>]   = dQclengthscale<EOL>dPinf[:,:,<NUM_LIT:0>] = dPinf_variance<EOL>dPinf[:,:,<NUM_LIT:1>] = dPinf_lengthscale<EOL>P0 = Pinf.copy()<EOL>dP0 = dPinf.copy()<EOL>if self.balance:<EOL><INDENT>import GPy.models.state_space_main as ssm<EOL>(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 )<EOL><DEDENT>return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)<EOL>
Return the state space representation of the covariance.

Note! For Sparse GP inference too small or two high values of lengthscale
lead to instabilities. This is because Qc are too high or too low
and P_inf are not full rank. This effect depends on approximatio order.
For N = 10. lengthscale must be in (0.8,8). For other N tests must be conducted.
N=6: (0.06,31)
Variance should be within reasonable bounds as well, but its dependence is linear.

The above facts do not take into accout regularization.