diff --git a/NEWS b/NEWS
index 32653cb..356439f 100644
--- a/NEWS
+++ b/NEWS
@@ -22,7 +22,11 @@ The following classes have been renamed:
 - none
 
 The following methods have been added:
-- none
+- GVector GModel::eval(GObservation&, GMatrixSparse*)
+- GVector GModelSky::eval(GObservation&, GMatrixSparse*)
+- GVector GResponse::convolve(GModelSky&, GObservation&, GMatrixSparse*)
+- GVector GObservation::model(GModels&, GMatrixSparse*)
+- GVector GObservation::model_grad(GModel&, GModelPar&) 
 
 The following methods have been removed:
 - none
diff --git a/include/GModel.hpp b/include/GModel.hpp
index 4b73c42..7979069 100644
--- a/include/GModel.hpp
+++ b/include/GModel.hpp
@@ -37,6 +37,8 @@
 #include "GTime.hpp"
 
 /* __ Forward declarations _______________________________________________ */
+class GVector;
+class GMatrixSparse;
 class GEvent;
 class GObservation;
 
@@ -125,6 +127,10 @@ public:
     virtual void        write(GXmlElement& xml) const = 0;
     virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
 
+    // Virtual methods
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse* gradients = NULL) const;
+
     // Implemented methods
     int                     size(void) const;
     int                     scales(void) const;
diff --git a/include/GModelData.hpp b/include/GModelData.hpp
index 8344f14..7b29d04 100644
--- a/include/GModelData.hpp
+++ b/include/GModelData.hpp
@@ -77,7 +77,7 @@ public:
     virtual void        read(const GXmlElement& xml) = 0;
     virtual void        write(GXmlElement& xml) const = 0;
     virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
-    
+
 protected:
     // Protected methods
     void init_members(void);
diff --git a/include/GModelSky.hpp b/include/GModelSky.hpp
index 4c35f96..55c9f20 100644
--- a/include/GModelSky.hpp
+++ b/include/GModelSky.hpp
@@ -1,7 +1,7 @@
 /***************************************************************************
  *                     GModelSky.hpp - Sky model class                     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -156,6 +156,10 @@ public:
     virtual void        write(GXmlElement& xml) const;
     virtual std::string print(const GChatter& chatter = NORMAL) const;
 
+    // Reimplemented virtual base class methods
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse*      gradients = NULL) const;
+
     // Other methods
     GModelSpatial*      spatial(void) const;
     GModelSpectral*     spectral(void) const;
diff --git a/include/GObservation.hpp b/include/GObservation.hpp
index 457b639..3f2c14d 100644
--- a/include/GObservation.hpp
+++ b/include/GObservation.hpp
@@ -102,6 +102,8 @@ public:
     virtual double           model(const GModels& models,
                                    const GEvent&  event,
                                    GVector*       gradient = NULL) const;
+    virtual GVector          model(const GModels& models,
+                                   GMatrixSparse* gradient = NULL) const;
     virtual int              nobserved(void) const;
     virtual double           npred(const GModels& models,
                                    GVector*       gradient = NULL) const;
@@ -109,14 +111,12 @@ public:
     virtual double           model_grad(const GModel&    model,
                                         const GModelPar& par,
                                         const GEvent&    event) const;
+    virtual GVector          model_grad(const GModel&    model,
+                                        const GModelPar& par) const;
     virtual double           npred_grad(const GModel&    model,
                                         const GModelPar& par) const;
     virtual void             remove_response_cache(const std::string& name);
 
-    // New methods
-    virtual GVector          model(const GModels& models,
-                                   GMatrixSparse* gradient = NULL) const;
-
     // Implemented methods
     bool               has_events(void) const;
     void               name(const std::string& name);
diff --git a/include/GResponse.hpp b/include/GResponse.hpp
index 7f81fc7..16bf796 100644
--- a/include/GResponse.hpp
+++ b/include/GResponse.hpp
@@ -37,7 +37,9 @@
 #include "GNdarray.hpp"
 
 /* __ Forward declarations _______________________________________________ */
+class GVector;
 class GMatrix;
+class GMatrixSparse;
 class GEvent;
 class GPhoton;
 class GSource;
@@ -104,6 +106,9 @@ public:
                                  const GEvent&       event,
                                  const GObservation& obs,
                                  const bool&         grad = true) const;
+    virtual GVector     convolve(const GModelSky&    model,
+                                 const GObservation& obs,
+                                 GMatrixSparse*      gradients = NULL) const;
     virtual double      irf_spatial(const GEvent&       event,
                                     const GSource&      source,
                                     const GObservation& obs) const;
diff --git a/src/model/GModel.cpp b/src/model/GModel.cpp
index ae88b3d..029eb8e 100644
--- a/src/model/GModel.cpp
+++ b/src/model/GModel.cpp
@@ -31,9 +31,13 @@
 #include "GTools.hpp"
 #include "GException.hpp"
 #include "GModel.hpp"
+#include "GVector.hpp"
+#include "GMatrixSparse.hpp"
+#include "GObservation.hpp"
 
 /* __ Method name definitions ____________________________________________ */
 #define G_ACCESS                           "GModel::operator[](std::string&)"
+#define G_EVAL                  "GModel::eval(GObservation&, GMatrixSparse*)"
 #define G_AT                                               "GModel::at(int&)"
 #define G_SCALE                              "GModelPar& GModel::scale(int&)"
 #define G_WRITE_SCALES                   "GModel::write_scales(GXmlElement&)"
@@ -221,6 +225,91 @@ const GModelPar& GModel::operator[](const std::string& name) const
  ==========================================================================*/
 
 /***********************************************************************//**
+ * @brief Return model values and gradients
+ *
+ * @param[in] obs Observation.
+ * @param[out] gradients Pointer to matrix of gradients.
+ * @return Model values.
+ *
+ * Evaluates the model values and parameter gradients for all events in an
+ * observation. Gradients are only returned if the @p gradients pointer is
+ * not NULL.
+ *
+ * The matrix of gradients is a sparse matrix where the number of rows
+ * corresponds to the number of model parameters and the number of columns
+ * corresponds to the number of events. An exception is thrown if the
+ * dimension of the @p gradients matrix is not compatible with the model
+ * and the observations.
+ ***************************************************************************/
+GVector GModel::eval(const GObservation& obs,
+                     GMatrixSparse*      gradients) const
+{
+    // Get number of model parameters and number of events
+    int npars   = size();
+    int nevents = obs.events()->size();
+
+    // Initialise gradients flag
+    bool grad = (gradients != NULL);
+
+    // Check matrix consistency
+    if (grad) {
+        if (gradients->rows() != npars) {
+            std::string msg = "Number of "+gammalib::str(gradients->rows())+
+                              " rows in gradient matrix differs from number "
+                              "of "+gammalib::str(npars)+" parameters "
+                              "in model. Please provide a compatible gradient "
+                              "matrix.";
+            throw GException::invalid_argument(G_EVAL, msg);
+        }
+        if (gradients->columns() != nevents) {
+            std::string msg = "Number of "+gammalib::str(gradients->columns())+
+                              " columns in gradient matrix differs from number "
+                              "of "+gammalib::str(nevents)+" events in "
+                              "observation. Please provide a compatible "
+                              "gradient matrix.";
+            throw GException::invalid_argument(G_EVAL, msg);
+        }
+    }
+
+    // Initialise values
+    GVector values(nevents);
+
+    // Loop over events
+    for (int k = 0; k < nevents; ++k) {
+
+        // Get reference to event
+        const GEvent& event = *((*obs.events())[k]);
+
+        // Evaluate probability
+        values[k] = eval(event, obs, grad);
+
+        // Set gradients
+        if (grad) {
+
+            // Initialise gradient vector
+            GVector gradient(npars);
+
+            // Extract relevant parameter gradients
+            for (int i = 0; i < npars; ++i) {
+                const GModelPar& par = (*this)[i];
+                if (par.is_free() && par.has_grad()) {
+                    gradient[i] = par.factor_gradient();
+                }
+            }
+
+            // Insert gradient vector into matrix
+            gradients->add_to_column(k, gradient);
+
+        } // endif: looped over gradients
+
+    } // endfor: looped over events
+
+    // Return values
+    return values;
+}
+
+
+/***********************************************************************//**
  * @brief Returns reference to model parameter by index
  *
  * @param[in] index Parameter index [0,...,size()-1].
@@ -298,6 +387,8 @@ bool GModel::has_par(const std::string& name) const
 /***********************************************************************//**
  * @brief Returns instruments to which model applies
  *
+ * @return Instruments.
+ *
  * Returns a comma separated list of instruments to which model applies. If
  * no instrument exists then an empty string is returned.
  ***************************************************************************/
@@ -485,6 +576,8 @@ void GModel::scale(const GModelPar& par)
 /***********************************************************************//**
  * @brief Returns observation identifiers to which model applies
  *
+ * @return Observation identifiers.
+ *
  * Returns a comma separated list of observation identifiers to which model
  * applies. If no observation identifier exists then an empty string is
  * returned.
@@ -542,7 +635,7 @@ void GModel::ids(const std::string& ids)
  *
  * @param[in] instrument Instrument name.
  * @param[in] id Observation identifier.
- * @return Validity flag
+ * @return Validity flag.
  *
  * Checks if specified instrument name and observation identifier is in list
  * of applicable instruments and identifiers. The check is case sensitive.
diff --git a/src/model/GModelSky.cpp b/src/model/GModelSky.cpp
index c2db277..db4d842 100644
--- a/src/model/GModelSky.cpp
+++ b/src/model/GModelSky.cpp
@@ -574,6 +574,30 @@ double GModelSky::eval(const GEvent&       event,
 
 
 /***********************************************************************//**
+ * @brief Evaluate sky model for all events of an observation
+ *
+ * @param[in] obs Observation.
+ * @param[out] gradients Pointer to sparse matrix holding model gradients.
+ * @return Values of sky model.
+ *
+ * Evalutes the value of the sky model for an @p event of a specific
+ * observation @p obs.
+ *
+ * If the @p gradients flag is true the method will also compute the
+ * parameter gradients for all model parameters.
+ ***************************************************************************/
+GVector GModelSky::eval(const GObservation& obs,
+                        GMatrixSparse*      gradients) const
+{
+    // Evaluate model
+    GVector values = obs.response()->convolve(*this, obs, gradients);
+
+    // Return
+    return values;
+}
+
+
+/***********************************************************************//**
  * @brief Return spatially integrated sky model
  *
  * @param[in] obsEng Measured photon energy.
@@ -595,7 +619,7 @@ double GModelSky::eval(const GEvent&       event,
  *                  S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  * \f]
  *
- * where                         
+ * where
  * \f$S(p,E,t)\f$ is the source model,
  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  * \f$p'\f$ is the measured photon direction,
@@ -1118,7 +1142,7 @@ void GModelSky::copy_members(const GModelSky& model)
 {
     // Copy attributes
     m_type = model.m_type;
-    
+
     // Clone model components
     m_spatial  = (model.m_spatial  != NULL) ? model.m_spatial->clone()  : NULL;
     m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
@@ -1218,12 +1242,12 @@ void GModelSky::set_type(void)
 
     // Continue only if we have a spatial model component
     if (m_spatial != NULL) {
-    
+
         // Is spatial model a point source?
         if (dynamic_cast<const GModelSpatialPointSource*>(m_spatial) != NULL) {
             m_type = "PointSource";
         }
-        
+
         // ... otherwise, is spatial model a radial source?
         else if (dynamic_cast<const GModelSpatialRadial*>(m_spatial) != NULL) {
             m_type = "ExtendedSource";
@@ -1243,7 +1267,7 @@ void GModelSky::set_type(void)
         else {
             m_type = "DiffuseSource";
         }
-    
+
     } // endif: there was a spatial model component
 
     // Return
diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp
index 72f4d1e..36cdaa9 100644
--- a/src/obs/GObservation.cpp
+++ b/src/obs/GObservation.cpp
@@ -48,8 +48,8 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_LIKELIHOOD           "GObservation::likelihood(GModels&, GVector*,"\
                                                   " GMatrixSparse*, double*)"
-#define G_MODEL                   "GObservation::model(GModels&, GPointing&,"\
-                                    " GInstDir&, GEnergy&, GTime&, GVector*)"
+#define G_MODEL1           "GObservation::model(GModels&, GEvent&, GVector*)"
+#define G_MODEL2              "GObservation::model(GModels&, GMatrixSparse*)"
 #define G_EVENTS                                     "GObservation::events()"
 #define G_NPRED                                "GObservation::npred(GModel&)"
 #define G_NPRED_SPEC              "GObservation::npred_spec(GModel&, GTime&)"
@@ -248,6 +248,7 @@ double GObservation::likelihood(const GModels& models,
  * @param[in] models Model container.
  * @param[in] event Observed event.
  * @param[out] gradient Pointer to gradient vector (optional).
+ * @return Model value.
  *
  * @exception GException::invalid_value
  *            Dimension of gradient vector mismatches number of parameters.
@@ -306,7 +307,7 @@ double GObservation::model(const GModels& models,
                             gammalib::str(models.npars())+" elements "
                             "requested while vector only contains "+
                             gammalib::str(gradient->size())+" elements.";
-                        throw GException::invalid_value(G_MODEL, msg);
+                        throw GException::invalid_value(G_MODEL1, msg);
                     }
                     #endif
 
@@ -390,41 +391,55 @@ double GObservation::model(const GModels& models,
  * @brief Return model values and (optionally) gradients
  *
  * @param[in] models Model container.
- * @param[out] gradient Pointer to gradient vector (optional).
+ * @param[out] gradient Pointer to sparse gradient matrix.
+ * @return Vector of model values.
  *
  * @exception GException::invalid_argument
  *            Gradient matrix mismatches number of events or model parameters.
+ *
+ * Returns the model values for each event in the observation.
+ *
+ * If @p gradient is not NULL, the matrix contains on output the model
+ * factor gradients for all events. Each row of the @p gradient matrix
+ * corresponds to one event, the columns correspond to the parameters
+ * of the @p models container.
  ***************************************************************************/
-#define G_MODEL2 "GObservation::model(GModels&, GMatrixSparse*)"
 GVector GObservation::model(const GModels& models,
                             GMatrixSparse* gradient) const
 {
-    // Initialise model values
-    GVector model(events()->size());
-
-    // Set gradient usage flag and initialise gradient counter
+    // Initialise variables
+    int  nevents  = events()->size();
     bool use_grad = (gradient != NULL);
-    int  igrad    = 0;      // Reset gradient counter
+    int  igrad    = 0;
 
-    // If gradient is available then gradient size
+    // Initialise model values
+    GVector values(nevents);
+
+    // If gradient is available then check gradient size and initialise
+    // sparse matrix stack
     if (use_grad) {
-        if (gradient->columns() != models.size()) {
+
+        // Check number of columns
+        if (gradient->columns() != models.npars()) {
             std::string msg = "Number of "+gammalib::str(gradient->columns())+
                               " columns in gradient matrix differs from number "
-                              "of "+gammalib::str(models.size())+" parameters "
-                              "in model container Please specify compatible "
+                              "of "+gammalib::str(models.npars())+" parameters "
+                              "in model container. Please specify compatible "
                               "arguments.";
             throw GException::invalid_argument(G_MODEL2, msg);
         }
-        if (gradient->rows() != events()->size()) {
+
+        // Check number of rows
+        if (gradient->rows() != nevents) {
             std::string msg = "Number of "+gammalib::str(gradient->rows())+
                               " rows in gradient matrix differs from number "
-                              "of "+gammalib::str(events()->size())+" events "
+                              "of "+gammalib::str(nevents)+" events "
                               "in observation. Please specify compatible "
                               "arguments.";
             throw GException::invalid_argument(G_MODEL2, msg);
         }
-    }
+
+    } // endif: Gradient requested
 
     // Loop over models
     for (int i = 0; i < models.size(); ++i) {
@@ -437,83 +452,91 @@ GVector GObservation::model(const GModels& models,
             // observation identifier
             if (mptr->is_valid(instrument(), id())) {
 
-                // Loop over events
-                for (int k = 0; k < events()->size(); ++k) {
+                // If no gradients should be used then evaluate model and
+                // add it to values
+                if (!use_grad) {
+                    values += mptr->eval(*this, NULL);
+                }
 
-                    // Get reference to event
-                    const GEvent& event = *((*events())[k]);
+                // ... otherwise evaluate model and gradients and add to
+                // model vector and gradient matrix
+                else {
 
-                    // Compute value and add to model
-                    model[k] += mptr->eval(event, *this, use_grad);
+                    // Allocate gradient matrix
+                    GMatrixSparse gradients(mptr->size(), nevents);
 
-                    // Optionally determine model gradients
-                    if (use_grad) {
+                    // Initialise sparse matrix stack
+                    gradients.stack_init(mptr->size()*1000, 1000);
 
-                        // If model provides the parameter indices that were
-                        // updated by the eval() method then use them as this
-                        // is less time consuming than looping over all indices
-                        if (mptr->has_eval_indices()) {
+                    // Evaluate model and add to values
+                    values += mptr->eval(*this, &gradients);
 
-                            // Get number of relevant parameters
-                            int npars = mptr->eval_indices().size();
+                    // Destroy sparse matrix stack (fills stack in matrix)
+                    gradients.stack_destroy();
 
-                            // Loop over all relevant parameters
-                            for (int index = 0; index < npars; ++index) {
+                    // If model provides the parameter indices that were
+                    // updated by the eval() method then use them as this
+                    // is less time consuming than looping over all indices
+                    if (mptr->has_eval_indices()) {
 
-                                // Retrieve parameter index
-                                int ipar = mptr->eval_indices()[index];
+                        // Get number of relevant parameters
+                        int npars = mptr->eval_indices().size();
 
-                                // Get reference to model parameter
-                                const GModelPar& par = (*mptr)[ipar];
+                        // Loop over all relevant parameters
+                        for (int index = 0; index < npars; ++index) {
+
+                            // Retrieve parameter index
+                            int ipar = mptr->eval_indices()[index];
 
-                                // Determine gradient
-                                double grad = 0.0;
-                                if (par.is_free()) {
-                                    if (par.has_grad()) {
-                                        grad = par.factor_gradient();
-                                    }
-                                    else {
-                                        grad = model_grad(*mptr, par, event);
-                                    }
+                            // Get reference to model parameter
+                            const GModelPar& par = (*mptr)[ipar];
+
+                            // Determine gradient
+                            GVector grad(nevents);
+                            if (par.is_free()) {
+                                if (par.has_grad()) {
+                                    grad = gradients.row(ipar);
+                                }
+                                else {
+                                    grad = model_grad(*mptr, par);
                                 }
+                            }
 
-                                // Set gradient
-                                (*gradient)(k,igrad+ipar) = grad;
+                            // Set gradient
+                            gradient->column(igrad+ipar, grad);
 
-                            } // endfor: looped over all parameter indices
+                        } // endfor: looped over all parameter indices
 
-                        } // endif: parameter indices were available
+                    } // endif: parameter indices were available
 
-                        // ... otherwise loop over all parameter indices
-                        else {
+                    // ... otherwise loop over all parameter indices
+                    else {
 
-                            // Loop over all parameters
-                            for (int ipar = 0; ipar < mptr->size(); ++ipar) {
+                        // Loop over all parameters
+                        for (int ipar = 0; ipar < mptr->size(); ++ipar) {
 
-                                // Get reference to model parameter
-                                const GModelPar& par = (*mptr)[ipar];
+                            // Get reference to model parameter
+                            const GModelPar& par = (*mptr)[ipar];
 
-                                // Determine gradient
-                                double grad = 0.0;
-                                if (par.is_free()) {
-                                    if (par.has_grad()) {
-                                        grad = par.factor_gradient();
-                                    }
-                                    else {
-                                        grad = model_grad(*mptr, par, event);
-                                    }
+                            // Determine gradient
+                            GVector grad(nevents);
+                            if (par.is_free()) {
+                                if (par.has_grad()) {
+                                    grad = gradients.row(ipar);
                                 }
+                                else {
+                                    grad = model_grad(*mptr, par);
+                                }
+                            }
 
-                                // Set gradient
-                                (*gradient)(k,igrad+ipar) = grad;
-
-                            } // endfor: looped over all parameter indices
+                            // Set gradient
+                            gradient->column(igrad+ipar, grad);
 
-                        } // endelse: no parameter indices were available
+                        } // endfor: looped over all parameter indices
 
-                    } // endif: use model gradients
+                    } // endelse: no parameter indices were available
 
-                } // endfor: loop over events
+                } // endif: use model gradients
 
             } // endif: model component was valid for instrument
 
@@ -524,8 +547,8 @@ GVector GObservation::model(const GModels& models,
 
     } // endfor: Looped over models
 
-    // Return
-    return model;
+    // Return values
+    return values;
 }
 
 
@@ -819,6 +842,123 @@ double GObservation::model_grad(const GModel&    model,
 
 
 /***********************************************************************//**
+ * @brief Returns parameter gradients of model for all events
+ *
+ * @param[in] model Model.
+ * @param[in] par Model parameter.
+ *
+ * This method uses a robust but simple difference method to estimate
+ * parameter gradients that have not been provided by the model. We use here
+ * a simple method as this method is likely used for spatial model parameters,
+ * and the spatial model may eventually be noisy due to numerical integration
+ * limits.
+ *
+ * The step size for the simple method has been fixed to 0.0002, which
+ * corresponds to about 1 arcsec for parameters that are given in degrees.
+ * The reasoning behind this value is that parameters that use numerical
+ * gradients are typically angles, such as for example the position, and
+ * we want to achieve arcsec precision with this method.
+ ***************************************************************************/
+GVector GObservation::model_grad(const GModel&    model,
+                                 const GModelPar& par) const
+{
+    // Initialise gradients
+    GVector gradients(events()->size());
+
+    // Compute gradient only if parameter is free
+    if (par.is_free()) {
+
+        // Get non-const model pointer
+        GModelPar* ptr = const_cast<GModelPar*>(&par);
+
+        // Save current model parameter
+        GModelPar current = par;
+
+        // Get actual parameter value
+        double x = par.factor_value();
+
+        // Set fixed step size for computation of derivative.
+        // By default, the step size is fixed to 0.0002.
+        const double step_size = 0.0002; // ~1 arcsec
+        double       h         = step_size;
+
+        // Re-adjust the step-size h in case that the initial step size is
+        // larger than the allowed parameter range
+        if (par.has_min() && par.has_max()) {
+            double par_h = par.factor_max() - par.factor_min();
+            if (par_h < h) {
+                h = par_h;
+            }
+        }
+
+        // Continue only if step size is positive
+        if (h > 0.0) {
+
+            // Setup derivative function
+            //GObservation::model_func function(this, &model, ptr, &event);
+            //GDerivative              derivative(&function);
+
+            // If we are too close to the minimum boundary use a right sided
+            // difference ...
+            if (par.has_min() && ((x-par.factor_min()) < h)) {
+                //grad = derivative.right_difference(x, h);
+                // Compute fs1
+                ptr->factor_value(x);
+                GVector fs1 = model.eval(*this, NULL);
+
+                // Compute fs2
+                ptr->factor_value(x+h);
+                GVector fs2 = model.eval(*this, NULL);
+
+                // Compute derivative
+                gradients = (fs1 - fs2) / h;
+
+            }
+
+            // ... otherwise if we are too close to the maximum boundary use
+            // a left sided difference ...
+            else if (par.has_max() && ((par.factor_max()-x) < h)) {
+                //grad = derivative.left_difference(x, h);
+                // Compute fs1
+                ptr->factor_value(x);
+                GVector fs1 = model.eval(*this, NULL);
+
+                // Compute fs2
+                ptr->factor_value(x-h);
+                GVector fs2 = model.eval(*this, NULL);
+
+                // Compute derivative
+                gradients = (fs1 - fs2) / h;
+            }
+
+            // ... otherwise use a symmetric difference
+            else {
+                //grad = derivative.difference(x, h);
+                // Compute fs1
+                ptr->factor_value(x+h);
+                GVector fs1 = model.eval(*this, NULL);
+
+                // Compute fs2
+                ptr->factor_value(x-h);
+                GVector fs2 = model.eval(*this, NULL);
+
+                // Compute derivative
+                gradients = 0.5 * (fs1 - fs2) / h;
+            }
+
+        } // endif: step size was positive
+
+        // Restore current model parameter
+        *ptr = current;
+
+    } // endif: model parameter was free
+
+    // Return gradients
+    return gradients;
+}
+
+
+/***********************************************************************//**
  * @brief Returns parameter gradient of Npred
  *
  * @param[in] model Gamma-ray source model.
@@ -1159,7 +1299,7 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models,
 
         // Update Poissonian statistic (excluding factorial term for faster
         // computation)
-        value -= log(model);
+        value -= std::log(model);
 
         // Skip bin now if there are no non-zero derivatives
         if (ndev < 1) {
@@ -1255,10 +1395,10 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
     // Allocate some working arrays
     int*          inx    = new int[npars];
     double*       values = new double[npars];
-    GMatrixSparse wrk_grad(nevents, npars);
+    GMatrixSparse wrk_matrix(nevents, npars);
 
     // Compute model and derivative
-    GVector model_vector = this->model(models, &wrk_grad);
+    GVector model_vector = this->model(models, &wrk_matrix);
 
     // Iterate over all bins
     for (int i = 0; i < events()->size(); ++i) {
@@ -1307,15 +1447,15 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
         // Update Npred
         *npred += model;
 
-        // Multiply gradient by bin size
-        //wrk_grad *= bin->size();
+        // Extract working gradient multiplied by bin size
+        GVector wrk_grad = wrk_matrix.row(i) * bin->size();
 
         // Create index array of non-zero derivatives and initialise working
         // array
         int ndev = 0;
         for (int k = 0; k < npars; ++k) {
             values[k] = 0.0;
-            if (wrk_grad(i,k) != 0.0 && !gammalib::is_infinite(wrk_grad(i,k))) {
+            if (wrk_grad[k] != 0.0 && !gammalib::is_infinite(wrk_grad[k])) {
                 inx[ndev] = k;
                 ndev++;
             }
@@ -1330,7 +1470,7 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
 
             // Update Poissonian statistic (excluding factorial term for
             // faster computation)
-            value -= data * log(model) - model;
+            value -= data * std::log(model) - model;
 
             // Skip bin now if there are no non-zero derivatives
             if (ndev < 1) {
@@ -1347,8 +1487,8 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
 
                 // Initialise computation
                 register int jpar    = inx[jdev];
-                double       g       = wrk_grad(i,jpar) * bin->size();
-                double       fa_i    = fa * g * bin->size();
+                double       g       = wrk_grad[jpar];
+                double       fa_i    = fa * g;
 
                 // Update gradient
                 (*gradient)[jpar] += fc * g;
@@ -1356,7 +1496,7 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
                 // Loop over rows
                 register int* ipar = inx;
                 for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
-                    values[idev] = fa_i * wrk_grad(i,*ipar);
+                    values[idev] = fa_i * wrk_grad[*ipar];
                 }
 
                 // Add column to matrix
@@ -1386,7 +1526,7 @@ double GObservation::likelihood_poisson_binned(const GModels& models,
             // Update gradient
             register int* ipar = inx;
             for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
-                (*gradient)[*ipar] += wrk_grad(i,*ipar) * bin->size();
+                (*gradient)[*ipar] += wrk_grad[*ipar];
             }
 
         } // endif: data was 0
diff --git a/src/obs/GObservations.cpp b/src/obs/GObservations.cpp
index ee9cf4e..30bc163 100644
--- a/src/obs/GObservations.cpp
+++ b/src/obs/GObservations.cpp
@@ -1,7 +1,7 @@
 /***************************************************************************
  *                GObservations.cpp - Observation container class          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2019 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -31,6 +31,8 @@
 #include "GTools.hpp"
 #include "GException.hpp"
 #include "GFilename.hpp"
+#include "GXml.hpp"
+#include "GModel.hpp"
 #include "GObservations.hpp"
 #include "GObservationRegistry.hpp"
 #include "GMatrixSparse.hpp"
diff --git a/src/obs/GObservations_likelihood.cpp b/src/obs/GObservations_likelihood.cpp
index ee4f9b3..494e428 100644
--- a/src/obs/GObservations_likelihood.cpp
+++ b/src/obs/GObservations_likelihood.cpp
@@ -1,7 +1,7 @@
 /***************************************************************************
  *         GObservations_likelihood.cpp - Likelihood function class        *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -31,6 +31,7 @@
 #include "GObservations.hpp"
 #include "GTools.hpp"
 #include "GFilename.hpp"
+#include "GModel.hpp"
 #include "GEvent.hpp"
 #include "GEventList.hpp"
 #include "GEventCube.hpp"
diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp
index 58ee3bb..9f46233 100644
--- a/src/obs/GResponse.cpp
+++ b/src/obs/GResponse.cpp
@@ -35,6 +35,7 @@
 #include "GMath.hpp"
 #include "GVector.hpp"
 #include "GMatrix.hpp"
+#include "GMatrixSparse.hpp"
 #include "GIntegral.hpp"
 #include "GIntegrals.hpp"
 #include "GResponse.hpp"
@@ -56,6 +57,8 @@
 #include "GModelSpatialComposite.hpp"
 
 /* __ Method name definitions ____________________________________________ */
+#define G_CONVOLVE          "GResponse::convolve(GModelSky&, GObservation&, "\
+                                                            "GMatrixSparse*)"
 #define G_IRF_RADIAL               "GResponse::irf_radial(GEvent&, GSource&,"\
                                                             " GObservation&)"
 #define G_IRF_ELLIPTICAL       "GResponse::irf_elliptical(GEvent&, GSource&,"\
@@ -311,6 +314,254 @@ double GResponse::convolve(const GModelSky&    model,
 
 
 /***********************************************************************//**
+ * @brief Convolve sky model with the instrument response
+ *
+ * @param[in] model Sky model.
+ * @param[in] obs Observation.
+ * @param[out] gradients Pointer to matrix of gradients.
+ * @return Event probabilities.
+ *
+ * Computes the event probability
+ *
+ * \f[
+ *    P(p',E',t') = \int \int \int
+ *                  S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
+ * \f]
+ *
+ * without taking into account any time dispersion. Energy dispersion is
+ * correctly handled by this method. If time dispersion is indeed needed,
+ * an instrument specific method needs to be provided.
+ ***************************************************************************/
+GVector GResponse::convolve(const GModelSky&    model,
+                            const GObservation& obs,
+                            GMatrixSparse*      gradients) const
+{
+    // Set number of iterations for Romberg integration.
+    static const int iter = 6;
+
+    // Get number of model parameters and number of events
+    int npars   = model.size();
+    int nevents = obs.events()->size();
+
+    // Initialise gradients flag
+    bool grad = (gradients != NULL);
+
+    // Check matrix consistency
+    if (grad) {
+        if (gradients->rows() != npars) {
+            std::string msg = "Number of "+gammalib::str(gradients->rows())+
+                              " rows in gradient matrix differs from number "
+                              "of "+gammalib::str(npars)+" parameters "
+                              "in model. Please provide a compatible gradient "
+                              "matrix.";
+            throw GException::invalid_argument(G_CONVOLVE, msg);
+        }
+        if (gradients->columns() != nevents) {
+            std::string msg = "Number of "+gammalib::str(gradients->columns())+
+                              " columns in gradient matrix differs from number "
+                              "of "+gammalib::str(nevents)+" events in "
+                              "observation. Please provide a compatible "
+                              "gradient matrix.";
+            throw GException::invalid_argument(G_CONVOLVE, msg);
+        }
+    }
+
+    // Initialise result
+    GVector probs(nevents);
+
+    // Continue only if the model has a spatial component
+    if (model.spatial() != NULL) {
+
+        // Loop over events
+        for (int k = 0; k < nevents; ++k) {
+
+            // Get reference to event
+            const GEvent& event = *((*obs.events())[k]);
+
+            // Get source time (no dispersion)
+            GTime srcTime = event.time();
+
+            // Case A: Integration
+            if (use_edisp()) {
+
+                // Retrieve true energy boundaries
+                GEbounds ebounds = this->ebounds(event.energy());
+
+                // Initialise Ndarray array
+                GNdarray array;
+
+                // Loop over all boundaries
+                for (int i = 0; i < ebounds.size(); ++i) {
+
+                    // Get true energy boundaries in MeV
+                    double etrue_min = ebounds.emin(i).MeV();
+                    double etrue_max = ebounds.emax(i).MeV();
+
+                    // Continue only if valid
+                    if (etrue_max > etrue_min) {
+
+                        // Setup integration function
+                        edisp_kern integrand(this, &obs, &model, &event, srcTime, grad);
+                        GIntegrals integral(&integrand);
+
+                        // Set number of iterations
+                        integral.fixed_iter(iter);
+
+                        // Do Romberg integration
+                        if (array.size() == 0) {
+                            array = integral.romberg(etrue_min, etrue_max);
+                        }
+                        else {
+                            array += integral.romberg(etrue_min, etrue_max);
+                        }
+
+                    } // endif: interval was valid
+
+                } // endfor: looped over intervals
+
+                // Initialise array index
+                int index = 0;
+
+                // Get probability
+                probs[k] = array(index++);
+
+                // Set gradients
+                if (grad) {
+
+                    // Initialise gradient vector
+                    GVector gradient(npars);
+
+                    // Get number of spatial, spectral and temporal parameters
+                    int n_spat = model.spatial()->size();
+                    int n_spec = model.spectral()->size();
+                    int n_temp = model.temporal()->size();
+
+                    // Set spectral gradients
+                    if (model.spectral() != NULL) {
+                        int offset = n_spat;
+                        for (int i = 0; i < n_spec; ++i) {
+                            GModelPar& par = (*(model.spectral()))[i];
+                            if (par.is_free() && par.has_grad()) {
+                                gradient[offset+i] = array(index++);
+                            }
+                        }
+                    }
+
+                    // Set temporal gradients
+                    if (model.temporal() != NULL) {
+                        int offset = n_spat + n_spec;
+                        for (int i = 0; i < n_temp; ++i) {
+                            GModelPar& par = (*(model.temporal()))[i];
+                            if (par.is_free() && par.has_grad()) {
+                                gradient[offset+i] = array(index++);
+                            }
+                        }
+                    }
+
+                    // Optionally set scale gradient for instrument
+                    if (model.has_scales()) {
+                        int offset = n_spat + n_spec + n_temp;
+                        for (int i = 0; i < model.scales(); ++i) {
+                            const GModelPar& par = model.scale(i);
+                            if (par.name() == obs.instrument()) {
+                                if (par.is_free() && par.has_grad()) {
+                                    gradient[offset+i] = array(index++);
+                                }
+                            }
+                        }
+                    }
+
+                    // Insert gradient vector into matrix
+                    gradients->add_to_column(k, gradient);
+
+                } // endif: set gradients
+
+            } // endif: Case A
+
+            // Case B: No integration (assume no energy dispersion)
+            else {
+
+                // Get source energy (no dispersion)
+                GEnergy srcEng  = event.energy();
+
+                // Evaluate probability
+                probs[k] = eval_prob(model, event, srcEng, srcTime, obs, grad);
+
+                // Set gradients
+                if (grad) {
+
+                    // Initialise gradient vector
+                    GVector gradient(npars);
+
+                    // Get number of spatial, spectral and temporal parameters
+                    int n_spat = model.spatial()->size();
+                    int n_spec = model.spectral()->size();
+                    int n_temp = model.temporal()->size();
+
+                    // Set spectral gradients
+                    if (model.spectral() != NULL) {
+                        int offset = n_spat;
+                        for (int i = 0; i < n_spec; ++i) {
+                            GModelPar& par = (*(model.spectral()))[i];
+                            if (par.is_free() && par.has_grad()) {
+                                gradient[offset+i] = par.factor_gradient();
+                            }
+                        }
+                    }
+
+                    // Set temporal gradients
+                    if (model.temporal() != NULL) {
+                        int offset = n_spat + n_spec;
+                        for (int i = 0; i < n_temp; ++i) {
+                            GModelPar& par = (*(model.temporal()))[i];
+                            if (par.is_free() && par.has_grad()) {
+                                gradient[offset+i] = par.factor_gradient();
+                            }
+                        }
+                    }
+
+                    // Optionally set scale gradient for instrument
+                    if (model.has_scales()) {
+                        int offset = n_spat + n_spec + n_temp;
+                        for (int i = 0; i < model.scales(); ++i) {
+                            const GModelPar& par = model.scale(i);
+                            if (par.name() == obs.instrument()) {
+                                if (par.is_free() && par.has_grad()) {
+                                    gradient[offset+i] = par.factor_gradient();
+                                }
+                            }
+                        }
+                    }
+
+                    // Insert gradient vector into matrix
+                    gradients->add_to_column(k, gradient);
+
+                } // endif: set gradients
+
+            } // endelse: Case B
+
+            // Compile option: Check for NaN/Inf
+            #if defined(G_NAN_CHECK)
+            if (gammalib::is_notanumber(probs[k]) || gammalib::is_infinite(probs[k])) {
+                std::cout << "*** ERROR: GResponse::convolve:";
+                std::cout << " NaN/Inf encountered";
+                std::cout << " (probs[" << k << "]=" << probs[k];
+                std::cout << ", event=" << event;
+                std::cout << ", srcTime=" << srcTime;
+                std::cout << ")" << std::endl;
+            }
+            #endif
+
+        } // endfor: looped over events
+
+    } // endif: spatial component valid
+
+    // Return probabilities
+    return probs;
+}
+
+
+/***********************************************************************//**
  * @brief Return instrument response integrated over the spatial model
  *
  * @param[in] event Event.