def sde(self): 
        """ 
        Return the state space representation of the standard periodic covariance.
        
        
        ! Note: one must constrain lengthscale not to drop below 0.2. (independently of approximation order)
        After this Bessel functions of the first becomes NaN. Rescaling
        time variable might help.
        
        ! Note: one must keep period also not very low. Because then
        the gradients wrt wavelength become ustable. 
        However this might depend on the data. For test example with
        300 data points the low limit is 0.15. 
        """ 
        
        #import pdb; pdb.set_trace()
        # Params to use: (in that order)
        #self.variance
        #self.period
        #self.lengthscale
        if self.approx_order is not None:
            N = int(self.approx_order)
        else:
            N = 7 # approximation order        
        
        p_period = float(self.period)        
        p_lengthscale = 2*float(self.lengthscale)
        p_variance = float(self.variance)        
        
        w0 = 2*np.pi/p_period # frequency
        # lengthscale is multiplied by 2 because of different definition of lengthscale
        
        [q2,dq2l] = seriescoeff(N, p_lengthscale, p_variance)        
        
        dq2l = 2*dq2l  # This is because the lengthscale if multiplied by 2.
        
        eps = 1e-12
        if np.any( np.isfinite(q2) == False) or np.any( np.abs(q2) > 1.0/eps) or np.any( np.abs(q2) < eps):
            warnings.warn("sde_Periodic:  Infinite, too small, or too large (eps={0:e}) values in q2 :".format(eps) + q2.__format__("") )
                                
        if np.any( np.isfinite(dq2l) == False) or np.any( np.abs(dq2l) > 1.0/eps) or np.any( np.abs(dq2l) < eps):
            warnings.warn("sde_Periodic:  Infinite, too small, or too large (eps={0:e}) values in dq2l :".format(eps) + q2.__format__("") )
                 
                 
        F    = np.kron(np.diag(range(0,N+1)),np.array( ((0, -w0), (w0, 0)) ) )
        L    = np.eye(2*(N+1))
        Qc   = np.zeros((2*(N+1), 2*(N+1)))
        P_inf = np.kron(np.diag(q2),np.eye(2))
        H    = np.kron(np.ones((1,N+1)),np.array((1,0)) )
        P0 = P_inf.copy()
        
        # Derivatives
        dF = np.empty((F.shape[0], F.shape[1], 3))
        dQc = np.empty((Qc.shape[0], Qc.shape[1], 3))
        dP_inf = np.empty((P_inf.shape[0], P_inf.shape[1], 3))         
        
        # Derivatives wrt self.variance
        dF[:,:,0] = np.zeros(F.shape)
        dQc[:,:,0] = np.zeros(Qc.shape)
        dP_inf[:,:,0] = P_inf / p_variance

        # Derivatives self.period
        dF[:,:,1] = np.kron(np.diag(range(0,N+1)),np.array( ((0,  w0), (-w0, 0)) ) / p_period );
        dQc[:,:,1] = np.zeros(Qc.shape)
        dP_inf[:,:,1] = np.zeros(P_inf.shape)      
        
        # Derivatives self.lengthscales        
        dF[:,:,2] = np.zeros(F.shape)
        dQc[:,:,2] = np.zeros(Qc.shape)
        dP_inf[:,:,2] = np.kron(np.diag(dq2l),np.eye(2))
        dP0 = dP_inf.copy()
        
        if self.balance:
            # Benefits of this are not very sound.
            import GPy.models.state_space_main as ssm
            (F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf,dP0) = ssm.balance_ss_model(F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0 )
            
        return (F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0)