Desctrete prediction function

        Input:
            k:int
                Iteration No. Starts at 0. Total number of iterations equal to the
                number of measurements.

            p_m: matrix of size (state_dim, time_series_no)
                Mean value from the previous step. For "multiple time series mode"
                it is matrix, second dimension of which correspond to different
                time series.

            p_P:
                Covariance matrix from the previous step.

            p_dyn_model_callable: class


            calc_grad_log_likelihood: boolean
                Whether to calculate gradient of the marginal likelihood
                of the state-space model. If true then the next parameter must
                provide the extra parameters for gradient calculation.

            p_dm: 3D array (state_dim, time_series_no, parameters_no)
                Mean derivatives from the previous step. For "multiple time series mode"
                it is 3D array, second dimension of which correspond to different
                time series.

            p_dP: 3D array (state_dim, state_dim, parameters_no)
                Mean derivatives from the previous step

        Output:
        ----------------------------
        m_pred, P_pred, dm_pred, dP_pred: metrices, 3D objects
            Results of the prediction steps.

def _kalman_prediction_step(k, p_m , p_P, p_dyn_model_callable, calc_grad_log_likelihood=False,
                                p_dm = None, p_dP = None):
        """
        Desctrete prediction function

        Input:
            k:int
                Iteration No. Starts at 0. Total number of iterations equal to the
                number of measurements.

            p_m: matrix of size (state_dim, time_series_no)
                Mean value from the previous step. For "multiple time series mode"
                it is matrix, second dimension of which correspond to different
                time series.

            p_P:
                Covariance matrix from the previous step.

            p_dyn_model_callable: class


            calc_grad_log_likelihood: boolean
                Whether to calculate gradient of the marginal likelihood
                of the state-space model. If true then the next parameter must
                provide the extra parameters for gradient calculation.

            p_dm: 3D array (state_dim, time_series_no, parameters_no)
                Mean derivatives from the previous step. For "multiple time series mode"
                it is 3D array, second dimension of which correspond to different
                time series.

            p_dP: 3D array (state_dim, state_dim, parameters_no)
                Mean derivatives from the previous step

        Output:
        ----------------------------
        m_pred, P_pred, dm_pred, dP_pred: metrices, 3D objects
            Results of the prediction steps.

        """

        # index correspond to values from previous iteration.
        A = p_dyn_model_callable.Ak(k,p_m,p_P) # state transition matrix (or Jacobian)
        Q = p_dyn_model_callable.Qk(k) # state noise matrix

        # Prediction step ->
        m_pred = p_dyn_model_callable.f_a(k, p_m, A) # predicted mean
        P_pred = A.dot(p_P).dot(A.T) + Q # predicted variance
        # Prediction step <-

        if calc_grad_log_likelihood:
            dA_all_params = p_dyn_model_callable.dAk(k) # derivatives of A wrt parameters
            dQ_all_params = p_dyn_model_callable.dQk(k) # derivatives of Q wrt parameters

            param_number = p_dP.shape[2]

            # p_dm, p_dP - derivatives form the previoius step
            dm_pred = np.empty(p_dm.shape)
            dP_pred = np.empty(p_dP.shape)

            for j in range(param_number):
                dA = dA_all_params[:,:,j]
                dQ = dQ_all_params[:,:,j]

                dP = p_dP[:,:,j]
                dm = p_dm[:,:,j]
                dm_pred[:,:,j] = np.dot(dA, p_m) + np.dot(A, dm)
                # prediction step derivatives for current parameter:

                dP_pred[:,:,j] = np.dot( dA ,np.dot(p_P, A.T))
                dP_pred[:,:,j] += dP_pred[:,:,j].T
                dP_pred[:,:,j] += np.dot( A ,np.dot(dP, A.T)) + dQ

                dP_pred[:,:,j] = 0.5*(dP_pred[:,:,j] + dP_pred[:,:,j].T) #symmetrize
        else:
            dm_pred = None
            dP_pred = None

        return m_pred, P_pred, dm_pred, dP_pred