Repository: DoubleML/doubleml-for-py (language: python)

Issue:
Added external predictions for PLR model
PLR extended by the option to use external predictions

### Description
The PLR model was extended by the possibility to calculate causal effects from externally generated predictions. This has been implemented along the lines of the IRM model.
As with the IRM model, the _nuisance_est method has been adapted so that it offers the possibility of passing external predictions in the fit method of the parent class.
In addition, the unit tests for comparing the values calculated manually and in the package have been extended. These now also compare the results between predictions generated with nuisance learner inside the package or predictions generated outside.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 41a409eb..d9ffb93c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,3 +29,4 @@ share/python-wheels/
 MANIFEST
 *.idea
 *.vscode
+.flake8
diff --git a/doubleml/_utils.py b/doubleml/_utils.py
index 345f914f..c2ed9d4e 100644
--- a/doubleml/_utils.py
+++ b/doubleml/_utils.py
@@ -333,3 +333,9 @@ def _var_est(psi, psi_deriv, apply_cross_fitting, smpls, is_cluster_data,
     sigma2_hat = np.multiply(scaling, gamma_hat)
 
     return sigma2_hat, var_scaling_factor
+
+
+def _cond_targets(target, cond_sample):
+    cond_target = target.astype(float)
+    cond_target[np.invert(cond_sample)] = np.nan
+    return cond_target
diff --git a/doubleml/double_ml.py b/doubleml/double_ml.py
index 082c6f3c..dd3547ad 100644
--- a/doubleml/double_ml.py
+++ b/doubleml/double_ml.py
@@ -64,6 +64,9 @@ def __init__(self,
         self._sensitivity_elements = None
         self._sensitivity_params = None
 
+        # initialize external predictions
+        self._external_predictions_implemented = False
+
         # check resampling specifications
         if not isinstance(n_folds, int):
             raise TypeError('The number of folds must be of int type. '
@@ -124,7 +127,7 @@ def __init__(self,
             self.draw_sample_splitting()
 
         # initialize arrays according to obj_dml_data and the resampling settings
-        self._psi, self._psi_deriv, self._psi_elements,\
+        self._psi, self._psi_deriv, self._psi_elements, \
             self._coef, self._se, self._all_coef, self._all_se, self._all_dml1_coef = self._initialize_arrays()
 
         # also initialize bootstrap arrays with the default number of bootstrap replications
@@ -486,7 +489,7 @@ def __psi_deriv(self):
     def __all_se(self):
         return self._all_se[self._i_treat, self._i_rep]
 
-    def fit(self, n_jobs_cv=None, store_predictions=True, store_models=False):
+    def fit(self, n_jobs_cv=None, store_predictions=True, external_predictions=None, store_models=False):
         """
         Estimate DoubleML models.
 
@@ -505,6 +508,13 @@ def fit(self, n_jobs_cv=None, store_predictions=True, store_models=False):
             to analyze the fitted models or extract information like variable importance.
             Default is ``False``.
 
+        external_predictions : None or dict
+            If `None` all models for the learners are fitted and evaluated. If a dictionary containing predictions
+            for a specific learner is supplied, the model will use the supplied nuisance predictions instead. Has to
+            be a nested dictionary where the keys refer to the treatment and the keys of the nested dictionarys refer to the
+            corresponding learners.
+            Default is `None`.
+
         Returns
         -------
         self : object
@@ -523,6 +533,13 @@ def fit(self, n_jobs_cv=None, store_predictions=True, store_models=False):
             raise TypeError('store_models must be True or False. '
                             f'Got {str(store_models)}.')
 
+        # check if external predictions are implemented
+        if self._external_predictions_implemented:
+            # check prediction format
+            self._check_external_predictions(external_predictions)
+        elif not self._external_predictions_implemented and external_predictions is not None:
+            raise NotImplementedError(f"External predictions not implemented for {self.__class__.__name__}.")
+
         # initialize rmse arrays for nuisance functions evaluation
         self._initialize_rmses()
 
@@ -546,8 +563,24 @@ def fit(self, n_jobs_cv=None, store_predictions=True, store_models=False):
                 if self._dml_data.n_treat > 1:
                     self._dml_data.set_x_d(self._dml_data.d_cols[i_d])
 
+                # set the supplied predictions for the treatment and each learner (including None)
+                ext_prediction_dict = {}
+                for learner in self.params_names:
+                    if external_predictions is None:
+                        ext_prediction_dict[learner] = None
+                    elif learner in external_predictions[self._dml_data.d_cols[i_d]].keys():
+                        if isinstance(external_predictions[self._dml_data.d_cols[i_d]][learner], np.ndarray):
+                            ext_prediction_dict[learner] = external_predictions[self._dml_data.d_cols[i_d]][learner][:, i_rep]
+                        else:
+                            ext_prediction_dict[learner] = None
+                    else:
+                        ext_prediction_dict[learner] = None
+
                 # ml estimation of nuisance models and computation of score elements
-                score_elements, preds = self._nuisance_est(self.__smpls, n_jobs_cv, return_models=store_models)
+                score_elements, preds = self._nuisance_est(self.__smpls, n_jobs_cv,
+                                                           external_predictions=ext_prediction_dict,
+                                                           return_models=store_models)
+
                 self._set_score_elements(score_elements, self._i_rep, self._i_treat)
 
                 # calculate rmses and store predictions and targets of the nuisance models
@@ -985,7 +1018,7 @@ def _initialize_ml_nuisance_params(self):
         pass
 
     @abstractmethod
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models):
+    def _nuisance_est(self, smpls, n_jobs_cv, return_models, external_predictions):
         pass
 
     @abstractmethod
@@ -1037,6 +1070,48 @@ def _check_learner(learner, learner_name, regressor, classifier):
 
         return learner_is_classifier
 
+    def _check_external_predictions(self, external_predictions):
+        if external_predictions is not None:
+            if not isinstance(external_predictions, dict):
+                raise TypeError('external_predictions must be a dictionary. '
+                                f'{str(external_predictions)} of type {str(type(external_predictions))} was passed.')
+
+            supplied_treatments = list(external_predictions.keys())
+            valid_treatments = self._dml_data.d_cols
+            if not set(supplied_treatments).issubset(valid_treatments):
+                raise ValueError('Invalid external_predictions. '
+                                 f'Invalid treatment variable in {str(supplied_treatments)}. '
+                                 'Valid treatment variables ' + ' or '.join(valid_treatments) + '.')
+
+            for treatment in supplied_treatments:
+                if not isinstance(external_predictions[treatment], dict):
+                    raise TypeError('external_predictions must be a nested dictionary. '
+                                    f'For treatment {str(treatment)} a value of type '
+                                    f'{str(type(external_predictions[treatment]))} was passed.')
+
+                supplied_learners = list(external_predictions[treatment].keys())
+                valid_learners = self.params_names
+                if not set(supplied_learners).issubset(valid_learners):
+                    raise ValueError('Invalid external_predictions. '
+                                     f'Invalid nuisance learner for treatment {str(treatment)} in {str(supplied_learners)}. '
+                                     'Valid nuisance learners ' + ' or '.join(valid_learners) + '.')
+
+                for learner in supplied_learners:
+                    if not isinstance(external_predictions[treatment][learner], np.ndarray):
+                        raise TypeError('Invalid external_predictions. '
+                                        'The values of the nested list must be a numpy array. '
+                                        'Invalid predictions for treatment ' + str(treatment) +
+                                        ' and learner ' + str(learner) + '. ' +
+                                        f'Object of type {str(type(external_predictions[treatment][learner]))} was passed.')
+
+                    expected_shape = (self._dml_data.n_obs, self.n_rep)
+                    if external_predictions[treatment][learner].shape != expected_shape:
+                        raise ValueError('Invalid external_predictions. '
+                                         f'The supplied predictions have to be of shape {str(expected_shape)}. '
+                                         'Invalid predictions for treatment ' + str(treatment) +
+                                         ' and learner ' + str(learner) + '. ' +
+                                         f'Predictions of shape {str(external_predictions[treatment][learner].shape)} passed.')
+
     def _initialize_arrays(self):
         # scores
         psi = np.full((self._dml_data.n_obs, self.n_rep, self._dml_data.n_coefs), np.nan)
diff --git a/doubleml/double_ml_cvar.py b/doubleml/double_ml_cvar.py
index f76c353c..7bc3fc7b 100644
--- a/doubleml/double_ml_cvar.py
+++ b/doubleml/double_ml_cvar.py
@@ -7,7 +7,7 @@
 from .double_ml import DoubleML
 from .double_ml_score_mixins import LinearScoreMixin
 from ._utils import _dml_cv_predict, _trimm, _predict_zero_one_propensity, \
-    _normalize_ipw, _dml_tune, _get_bracket_guess, _solve_ipw_score
+    _normalize_ipw, _dml_tune, _get_bracket_guess, _solve_ipw_score, _cond_targets
 from .double_ml_data import DoubleMLData
 from ._utils_resampling import DoubleMLResampling
 from ._utils_checks import _check_score, _check_trimming, _check_zero_one_treatment, _check_treatment, \
@@ -207,7 +207,7 @@ def _initialize_ml_nuisance_params(self):
         self._params = {learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols}
                         for learner in ['ml_g', 'ml_m']}
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
@@ -296,8 +296,7 @@ def ipw_score(theta):
         m_hat['targets'] = d
 
         # set the target for g to be a float and only relevant values
-        g_hat['targets'] = g_hat['targets'].astype(float)
-        g_hat['targets'][d != self.treatment] = np.nan
+        g_hat['targets'] = _cond_targets(g_hat['targets'], cond_sample=(d == self.treatment))
 
         if return_models:
             g_hat['models'] = fitted_models['ml_g']
diff --git a/doubleml/double_ml_did.py b/doubleml/double_ml_did.py
index 6d469bd3..77bface4 100644
--- a/doubleml/double_ml_did.py
+++ b/doubleml/double_ml_did.py
@@ -146,8 +146,8 @@ def __init__(self,
         self._trimming_rule = trimming_rule
         self._trimming_threshold = trimming_threshold
         _check_trimming(self._trimming_rule, self._trimming_threshold)
-
         self._sensitivity_implemented = True
+        self._external_predictions_implemented = True
 
     @property
     def in_sample_normalization(self):
@@ -194,7 +194,7 @@ def _check_data(self, obj_dml_data):
                              'needs to be specified as treatment variable.')
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
@@ -203,31 +203,49 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
         # nuisance g
         # get train indices for d == 0
         smpls_d0, smpls_d1 = _get_cond_smpls(smpls, d)
-        g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d0, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
-
-        _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat0['targets'] = g_hat0['targets'].astype(float)
-        g_hat0['targets'][d == 1] = np.nan
 
-        g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d1, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
+        # nuisance g for d==0
+        if external_predictions['ml_g0'] is not None:
+            g_hat0 = {'preds': external_predictions['ml_g0'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d0, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
+
+            _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            # adjust target values to consider only compatible subsamples
+            g_hat0['targets'] = g_hat0['targets'].astype(float)
+            g_hat0['targets'][d == 1] = np.nan
+
+        # nuisance g for d==1
+        if external_predictions['ml_g1'] is not None:
+            g_hat1 = {'preds': external_predictions['ml_g1'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d1, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
 
-        _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat1['targets'] = g_hat1['targets'].astype(float)
-        g_hat1['targets'][d == 0] = np.nan
+            _check_finite_predictions(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            # adjust target values to consider only compatible subsamples
+            g_hat1['targets'] = g_hat1['targets'].astype(float)
+            g_hat1['targets'][d == 0] = np.nan
 
         # only relevant for observational setting
         m_hat = {'preds': None, 'targets': None, 'models': None}
         if self.score == 'observational':
             # nuisance m
-            m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
-                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                    return_models=return_models)
+            if external_predictions['ml_m'] is not None:
+                m_hat = {'preds': external_predictions['ml_m'],
+                         'targets': None,
+                         'models': None}
+            else:
+                m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
+                                        est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                        return_models=return_models)
             _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
             _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
             m_hat['preds'] = _trimm(m_hat['preds'], self.trimming_rule, self.trimming_threshold)
diff --git a/doubleml/double_ml_did_cs.py b/doubleml/double_ml_did_cs.py
index 5591632b..55b5e32e 100644
--- a/doubleml/double_ml_did_cs.py
+++ b/doubleml/double_ml_did_cs.py
@@ -148,6 +148,7 @@ def __init__(self,
         _check_trimming(self._trimming_rule, self._trimming_threshold)
 
         self._sensitivity_implemented = True
+        self._external_predictions_implemented = True
 
     @property
     def in_sample_normalization(self):
@@ -207,7 +208,7 @@ def _check_data(self, obj_dml_data):
 
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
@@ -228,40 +229,62 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
 
         # nuisance g
         smpls_d0_t0, smpls_d0_t1, smpls_d1_t0, smpls_d1_t1 = _get_cond_smpls_2d(smpls, d, t)
-
-        g_hat_d0_t0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d0_t0, n_jobs=n_jobs_cv,
-                                      est_params=self._get_params('ml_g_d0_t0'), method=self._predict_method['ml_g'],
-                                      return_models=return_models)
-        g_hat_d0_t0['targets'] = g_hat_d0_t0['targets'].astype(float)
-        g_hat_d0_t0['targets'][np.invert((d == 0) & (t == 0))] = np.nan
-
-        g_hat_d0_t1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d0_t1, n_jobs=n_jobs_cv,
-                                      est_params=self._get_params('ml_g_d0_t1'), method=self._predict_method['ml_g'],
-                                      return_models=return_models)
-        g_hat_d0_t1['targets'] = g_hat_d0_t1['targets'].astype(float)
-        g_hat_d0_t1['targets'][np.invert((d == 0) & (t == 1))] = np.nan
-
-        g_hat_d1_t0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d1_t0, n_jobs=n_jobs_cv,
-                                      est_params=self._get_params('ml_g_d1_t0'), method=self._predict_method['ml_g'],
-                                      return_models=return_models)
-        g_hat_d1_t0['targets'] = g_hat_d1_t0['targets'].astype(float)
-        g_hat_d1_t0['targets'][np.invert((d == 1) & (t == 0))] = np.nan
-
-        g_hat_d1_t1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d1_t1, n_jobs=n_jobs_cv,
-                                      est_params=self._get_params('ml_g_d1_t1'), method=self._predict_method['ml_g'],
-                                      return_models=return_models)
-        g_hat_d1_t1['targets'] = g_hat_d1_t1['targets'].astype(float)
-        g_hat_d1_t1['targets'][np.invert((d == 1) & (t == 1))] = np.nan
+        if external_predictions['ml_g_d0_t0'] is not None:
+            g_hat_d0_t0 = {'preds': external_predictions['ml_g_d0_t0'],
+                           'targets': None,
+                           'models': None}
+        else:
+            g_hat_d0_t0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d0_t0, n_jobs=n_jobs_cv,
+                                          est_params=self._get_params('ml_g_d0_t0'), method=self._predict_method['ml_g'],
+                                          return_models=return_models)
+
+            g_hat_d0_t0['targets'] = g_hat_d0_t0['targets'].astype(float)
+            g_hat_d0_t0['targets'][np.invert((d == 0) & (t == 0))] = np.nan
+        if external_predictions['ml_g_d0_t1'] is not None:
+            g_hat_d0_t1 = {'preds': external_predictions['ml_g_d0_t1'],
+                           'targets': None,
+                           'models': None}
+        else:
+            g_hat_d0_t1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d0_t1, n_jobs=n_jobs_cv,
+                                          est_params=self._get_params('ml_g_d0_t1'), method=self._predict_method['ml_g'],
+                                          return_models=return_models)
+            g_hat_d0_t1['targets'] = g_hat_d0_t1['targets'].astype(float)
+            g_hat_d0_t1['targets'][np.invert((d == 0) & (t == 1))] = np.nan
+        if external_predictions['ml_g_d1_t0'] is not None:
+            g_hat_d1_t0 = {'preds': external_predictions['ml_g_d1_t0'],
+                           'targets': None,
+                           'models': None}
+        else:
+            g_hat_d1_t0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d1_t0, n_jobs=n_jobs_cv,
+                                          est_params=self._get_params('ml_g_d1_t0'), method=self._predict_method['ml_g'],
+                                          return_models=return_models)
+            g_hat_d1_t0['targets'] = g_hat_d1_t0['targets'].astype(float)
+            g_hat_d1_t0['targets'][np.invert((d == 1) & (t == 0))] = np.nan
+        if external_predictions['ml_g_d1_t1'] is not None:
+            g_hat_d1_t1 = {'preds': external_predictions['ml_g_d1_t1'],
+                           'targets': None,
+                           'models': None}
+        else:
+            g_hat_d1_t1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls_d1_t1, n_jobs=n_jobs_cv,
+                                          est_params=self._get_params('ml_g_d1_t1'), method=self._predict_method['ml_g'],
+                                          return_models=return_models)
+            g_hat_d1_t1['targets'] = g_hat_d1_t1['targets'].astype(float)
+            g_hat_d1_t1['targets'][np.invert((d == 1) & (t == 1))] = np.nan
 
         # only relevant for observational or experimental setting
         m_hat = {'preds': None, 'targets': None, 'models': None}
         if self.score == 'observational':
             # nuisance m
-            m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
-                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                    return_models=return_models)
-            _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
-            _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
+            if external_predictions['ml_m'] is not None:
+                m_hat = {'preds': external_predictions['ml_m'],
+                         'targets': None,
+                         'models': None}
+            else:
+                m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
+                                        est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                        return_models=return_models)
+                _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
+                _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
             m_hat['preds'] = _trimm(m_hat['preds'], self.trimming_rule, self.trimming_threshold)
 
         psi_a, psi_b = self._score_elements(y, d, t,
diff --git a/doubleml/double_ml_iivm.py b/doubleml/double_ml_iivm.py
index a872ebca..d981250e 100644
--- a/doubleml/double_ml_iivm.py
+++ b/doubleml/double_ml_iivm.py
@@ -193,6 +193,7 @@ def __init__(self,
                 raise TypeError("subgroups['never_takers'] must be True or False. "
                                 f'Got {str(subgroups["never_takers"])}.')
         self.subgroups = subgroups
+        self._external_predictions_implemented = True
 
     @property
     def normalize_ipw(self):
@@ -246,7 +247,7 @@ def _check_data(self, obj_dml_data):
             raise ValueError(err_msg)
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, z = check_X_y(x, np.ravel(self._dml_data.z),
@@ -258,13 +259,18 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
         smpls_z0, smpls_z1 = _get_cond_smpls(smpls, z)
 
         # nuisance g
-        g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_z0, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
-        _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat0['targets'] = g_hat0['targets'].astype(float)
-        g_hat0['targets'][z == 1] = np.nan
+        if external_predictions['ml_g0'] is not None:
+            g_hat0 = {'preds': external_predictions['ml_g0'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_z0, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
+            _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            # adjust target values to consider only compatible subsamples
+            g_hat0['targets'] = g_hat0['targets'].astype(float)
+            g_hat0['targets'][z == 1] = np.nan
 
         if self._dml_data.binary_outcome:
             binary_preds = (type_of_target(g_hat0['preds']) == 'binary')
@@ -276,14 +282,18 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
                                  'probabilities and not labels are predicted.')
 
             _check_is_propensity(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls, eps=1e-12)
-
-        g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_z1, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
-        _check_finite_predictions(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat1['targets'] = g_hat1['targets'].astype(float)
-        g_hat1['targets'][z == 0] = np.nan
+        if external_predictions['ml_g1'] is not None:
+            g_hat1 = {'preds': external_predictions['ml_g1'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_z1, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
+            _check_finite_predictions(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            # adjust target values to consider only compatible subsamples
+            g_hat1['targets'] = g_hat1['targets'].astype(float)
+            g_hat1['targets'][z == 0] = np.nan
 
         if self._dml_data.binary_outcome:
             binary_preds = (type_of_target(g_hat1['preds']) == 'binary')
@@ -297,34 +307,53 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
             _check_is_propensity(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls, eps=1e-12)
 
         # nuisance m
-        m_hat = _dml_cv_predict(self._learner['ml_m'], x, z, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                return_models=return_models)
-        _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
-        _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
+        if external_predictions['ml_m'] is not None:
+            m_hat = {'preds': external_predictions['ml_m'],
+                     'targets': None,
+                     'models': None}
+        else:
+            m_hat = _dml_cv_predict(self._learner['ml_m'], x, z, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                    return_models=return_models)
+            _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
+            _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
 
         # nuisance r
+        r0 = external_predictions['ml_r0'] is not None
         if self.subgroups['always_takers']:
-            r_hat0 = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls_z0, n_jobs=n_jobs_cv,
-                                     est_params=self._get_params('ml_r0'), method=self._predict_method['ml_r'],
-                                     return_models=return_models)
+            if r0:
+                r_hat0 = {'preds': external_predictions['ml_r0'],
+                          'targets': None,
+                          'models': None}
+            else:
+                r_hat0 = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls_z0, n_jobs=n_jobs_cv,
+                                         est_params=self._get_params('ml_r0'), method=self._predict_method['ml_r'],
+                                         return_models=return_models)
         else:
             r_hat0 = {'preds': np.zeros_like(d), 'targets': np.zeros_like(d), 'models': None}
-        _check_finite_predictions(r_hat0['preds'], self._learner['ml_r'], 'ml_r', smpls)
-        # adjust target values to consider only compatible subsamples
-        r_hat0['targets'] = r_hat0['targets'].astype(float)
-        r_hat0['targets'][z == 1] = np.nan
+        if not r0:
+            _check_finite_predictions(r_hat0['preds'], self._learner['ml_r'], 'ml_r', smpls)
+            # adjust target values to consider only compatible subsamples
+            r_hat0['targets'] = r_hat0['targets'].astype(float)
+            r_hat0['targets'][z == 1] = np.nan
 
+        r1 = external_predictions['ml_r1'] is not None
         if self.subgroups['never_takers']:
-            r_hat1 = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls_z1, n_jobs=n_jobs_cv,
-                                     est_params=self._get_params('ml_r1'), method=self._predict_method['ml_r'],
-                                     return_models=return_models)
+            if r1:
+                r_hat1 = {'preds': external_predictions['ml_r1'],
+                          'targets': None,
+                          'models': None}
+            else:
+                r_hat1 = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls_z1, n_jobs=n_jobs_cv,
+                                         est_params=self._get_params('ml_r1'), method=self._predict_method['ml_r'],
+                                         return_models=return_models)
         else:
             r_hat1 = {'preds': np.ones_like(d), 'targets': np.ones_like(d), 'models': None}
-        _check_finite_predictions(r_hat1['preds'], self._learner['ml_r'], 'ml_r', smpls)
-        # adjust target values to consider only compatible subsamples
-        r_hat1['targets'] = r_hat1['targets'].astype(float)
-        r_hat1['targets'][z == 0] = np.nan
+        if not r1:
+            _check_finite_predictions(r_hat1['preds'], self._learner['ml_r'], 'ml_r', smpls)
+            # adjust target values to consider only compatible subsamples
+            r_hat1['targets'] = r_hat1['targets'].astype(float)
+            r_hat1['targets'][z == 0] = np.nan
 
         psi_a, psi_b = self._score_elements(y, z, d,
                                             g_hat0['preds'], g_hat1['preds'], m_hat['preds'],
diff --git a/doubleml/double_ml_irm.py b/doubleml/double_ml_irm.py
index 4bbe42f7..9149462f 100644
--- a/doubleml/double_ml_irm.py
+++ b/doubleml/double_ml_irm.py
@@ -11,7 +11,7 @@
 from .double_ml_data import DoubleMLData
 from .double_ml_score_mixins import LinearScoreMixin
 
-from ._utils import _dml_cv_predict, _get_cond_smpls, _dml_tune, _trimm, _normalize_ipw
+from ._utils import _dml_cv_predict, _get_cond_smpls, _dml_tune, _trimm, _normalize_ipw, _cond_targets
 from ._utils_checks import _check_score, _check_trimming, _check_finite_predictions, _check_is_propensity, _check_integer
 
 
@@ -160,6 +160,7 @@ def __init__(self,
         _check_trimming(self._trimming_rule, self._trimming_threshold)
 
         self._sensitivity_implemented = True
+        self._external_predictions_implemented = True
 
     @property
     def normalize_ipw(self):
@@ -206,39 +207,51 @@ def _check_data(self, obj_dml_data):
                              'needs to be specified as treatment variable.')
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
                          force_all_finite=False)
         # get train indices for d == 0 and d == 1
         smpls_d0, smpls_d1 = _get_cond_smpls(smpls, d)
+        g0_external = external_predictions['ml_g0'] is not None
+        g1_external = external_predictions['ml_g1'] is not None
+        m_external = external_predictions['ml_m'] is not None
 
         # nuisance g
-        g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d0, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
-        _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat0['targets'] = g_hat0['targets'].astype(float)
-        g_hat0['targets'][d == 1] = np.nan
-
-        if self._dml_data.binary_outcome:
-            binary_preds = (type_of_target(g_hat0['preds']) == 'binary')
-            zero_one_preds = np.all((np.power(g_hat0['preds'], 2) - g_hat0['preds']) == 0)
-            if binary_preds & zero_one_preds:
-                raise ValueError(f'For the binary outcome variable {self._dml_data.y_col}, '
-                                 f'predictions obtained with the ml_g learner {str(self._learner["ml_g"])} are also '
-                                 'observed to be binary with values 0 and 1. Make sure that for classifiers '
-                                 'probabilities and not labels are predicted.')
-
-        g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d1, n_jobs=n_jobs_cv,
-                                 est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
-                                 return_models=return_models)
-        _check_finite_predictions(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls)
-        # adjust target values to consider only compatible subsamples
-        g_hat1['targets'] = g_hat1['targets'].astype(float)
-        g_hat1['targets'][d == 0] = np.nan
+        if g0_external:
+            # use external predictions
+            g_hat0 = {'preds': external_predictions['ml_g0'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat0 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d0, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g0'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
+            _check_finite_predictions(g_hat0['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            g_hat0['targets'] = _cond_targets(g_hat0['targets'], cond_sample=(d == 0))
+
+            if self._dml_data.binary_outcome:
+                binary_preds = (type_of_target(g_hat0['preds']) == 'binary')
+                zero_one_preds = np.all((np.power(g_hat0['preds'], 2) - g_hat0['preds']) == 0)
+                if binary_preds & zero_one_preds:
+                    raise ValueError(f'For the binary outcome variable {self._dml_data.y_col}, '
+                                     f'predictions obtained with the ml_g learner {str(self._learner["ml_g"])} are also '
+                                     'observed to be binary with values 0 and 1. Make sure that for classifiers '
+                                     'probabilities and not labels are predicted.')
+
+        if g1_external:
+            # use external predictions
+            g_hat1 = {'preds': external_predictions['ml_g1'],
+                      'targets': None,
+                      'models': None}
+        else:
+            g_hat1 = _dml_cv_predict(self._learner['ml_g'], x, y, smpls=smpls_d1, n_jobs=n_jobs_cv,
+                                     est_params=self._get_params('ml_g1'), method=self._predict_method['ml_g'],
+                                     return_models=return_models)
+            _check_finite_predictions(g_hat1['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            # adjust target values to consider only compatible subsamples
+            g_hat1['targets'] = _cond_targets(g_hat1['targets'], cond_sample=(d == 1))
 
         if self._dml_data.binary_outcome:
             binary_preds = (type_of_target(g_hat1['preds']) == 'binary')
@@ -250,12 +263,18 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
                                  'probabilities and not labels are predicted.')
 
         # nuisance m
-        m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                return_models=return_models)
-        _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
-        _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
-        m_hat['preds'] = _trimm(m_hat['preds'], self.trimming_rule, self.trimming_threshold)
+        if m_external:
+            # use external predictions
+            m_hat = {'preds': external_predictions['ml_m'],
+                     'targets': None,
+                     'models': None}
+        else:
+            m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                    return_models=return_models)
+            _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
+            _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
+            m_hat['preds'] = _trimm(m_hat['preds'], self.trimming_rule, self.trimming_threshold)
 
         psi_a, psi_b = self._score_elements(y, d,
                                             g_hat0['preds'], g_hat1['preds'], m_hat['preds'],
@@ -284,12 +303,15 @@ def _score_elements(self, y, d, g_hat0, g_hat1, m_hat, smpls):
             for _, test_index in smpls:
                 p_hat[test_index] = np.mean(d[test_index])
 
+        m_hat_adj = np.full_like(m_hat, np.nan, dtype='float64')
         if self.normalize_ipw:
             if self.dml_procedure == 'dml1':
                 for _, test_index in smpls:
-                    m_hat[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
+                    m_hat_adj[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
             else:
-                m_hat = _normalize_ipw(m_hat, d)
+                m_hat_adj = _normalize_ipw(m_hat, d)
+        else:
+            m_hat_adj = m_hat
 
         # compute residuals
         u_hat0 = y - g_hat0
@@ -300,19 +322,19 @@ def _score_elements(self, y, d, g_hat0, g_hat1, m_hat, smpls):
         if isinstance(self.score, str):
             if self.score == 'ATE':
                 psi_b = g_hat1 - g_hat0 \
-                    + np.divide(np.multiply(d, u_hat1), m_hat) \
-                    - np.divide(np.multiply(1.0-d, u_hat0), 1.0 - m_hat)
-                psi_a = np.full_like(m_hat, -1.0)
+                    + np.divide(np.multiply(d, u_hat1), m_hat_adj) \
+                    - np.divide(np.multiply(1.0-d, u_hat0), 1.0 - m_hat_adj)
+                psi_a = np.full_like(m_hat_adj, -1.0)
             else:
                 assert self.score == 'ATTE'
                 psi_b = np.divide(np.multiply(d, u_hat0), p_hat) \
-                    - np.divide(np.multiply(m_hat, np.multiply(1.0-d, u_hat0)),
-                                np.multiply(p_hat, (1.0 - m_hat)))
+                    - np.divide(np.multiply(m_hat_adj, np.multiply(1.0-d, u_hat0)),
+                                np.multiply(p_hat, (1.0 - m_hat_adj)))
                 psi_a = - np.divide(d, p_hat)
         else:
             assert callable(self.score)
             psi_a, psi_b = self.score(y=y, d=d,
-                                      g_hat0=g_hat0, g_hat1=g_hat1, m_hat=m_hat,
+                                      g_hat0=g_hat0, g_hat1=g_hat1, m_hat=m_hat_adj,
                                       smpls=smpls)
 
         return psi_a, psi_b
diff --git a/doubleml/double_ml_lpq.py b/doubleml/double_ml_lpq.py
index d4a75e82..2b7260eb 100644
--- a/doubleml/double_ml_lpq.py
+++ b/doubleml/double_ml_lpq.py
@@ -1,5 +1,4 @@
 import numpy as np
-import copy
 from sklearn.utils.multiclass import type_of_target
 from sklearn.base import clone
 from sklearn.utils import check_X_y
@@ -9,11 +8,19 @@
 from .double_ml_score_mixins import NonLinearScoreMixin
 from .double_ml_data import DoubleMLData
 
-from ._utils import _dml_cv_predict, _trimm, _predict_zero_one_propensity, \
-    _get_bracket_guess, _default_kde, _normalize_ipw, _dml_tune, _solve_ipw_score
+from ._utils import (
+    _dml_cv_predict,
+    _trimm,
+    _predict_zero_one_propensity,
+    _cond_targets,
+    _get_bracket_guess,
+    _default_kde,
+    _normalize_ipw,
+    _dml_tune,
+    _solve_ipw_score,
+)
 from ._utils_resampling import DoubleMLResampling
-from ._utils_checks import _check_score, _check_trimming, _check_zero_one_treatment, _check_treatment, \
-    _check_quantile
+from ._utils_checks import _check_score, _check_trimming, _check_zero_one_treatment, _check_treatment, _check_quantile
 
 
 class DoubleMLLPQ(NonLinearScoreMixin, DoubleML):
@@ -100,29 +107,25 @@ class DoubleMLLPQ(NonLinearScoreMixin, DoubleML):
     d  0.217244  0.636453  0.341336  0.73285 -1.03018  1.464668
     """
 
-    def __init__(self,
-                 obj_dml_data,
-                 ml_g,
-                 ml_m,
-                 treatment=1,
-                 quantile=0.5,
-                 n_folds=5,
-                 n_rep=1,
-                 score='LPQ',
-                 dml_procedure='dml2',
-                 normalize_ipw=True,
-                 kde=None,
-                 trimming_rule='truncate',
-                 trimming_threshold=1e-2,
-                 draw_sample_splitting=True,
-                 apply_cross_fitting=True):
-        super().__init__(obj_dml_data,
-                         n_folds,
-                         n_rep,
-                         score,
-                         dml_procedure,
-                         draw_sample_splitting,
-                         apply_cross_fitting)
+    def __init__(
+        self,
+        obj_dml_data,
+        ml_g,
+        ml_m,
+        treatment=1,
+        quantile=0.5,
+        n_folds=5,
+        n_rep=1,
+        score="LPQ",
+        dml_procedure="dml2",
+        normalize_ipw=True,
+        kde=None,
+        trimming_rule="truncate",
+        trimming_threshold=1e-2,
+        draw_sample_splitting=True,
+        apply_cross_fitting=True,
+    ):
+        super().__init__(obj_dml_data, n_folds, n_rep, score, dml_procedure, draw_sample_splitting, apply_cross_fitting)
 
         self._quantile = quantile
         self._treatment = treatment
@@ -130,21 +133,21 @@ def __init__(self,
             self._kde = _default_kde
         else:
             if not callable(kde):
-                raise TypeError('kde should be either a callable or None. '
-                                '%r was passed.' % kde)
+                raise TypeError("kde should be either a callable or None. " "%r was passed." % kde)
             self._kde = kde
         self._normalize_ipw = normalize_ipw
 
         self._check_data(self._dml_data)
 
-        valid_score = ['LPQ']
+        valid_score = ["LPQ"]
         _check_score(self.score, valid_score, allow_callable=False)
         _check_quantile(self.quantile)
         _check_treatment(self.treatment)
 
         if not isinstance(self.normalize_ipw, bool):
-            raise TypeError('Normalization indicator has to be boolean. ' +
-                            f'Object of type {str(type(self.normalize_ipw))} passed.')
+            raise TypeError(
+                "Normalization indicator has to be boolean. " + f"Object of type {str(type(self.normalize_ipw))} passed."
+            )
 
         # initialize starting values and bounds
         self._coef_bounds = (self._dml_data.y.min(), self._dml_data.y.max())
@@ -155,25 +158,36 @@ def __init__(self,
         self._trimming_threshold = trimming_threshold
         _check_trimming(self._trimming_rule, self._trimming_threshold)
 
-        _ = self._check_learner(ml_g, 'ml_g', regressor=False, classifier=True)
-        _ = self._check_learner(ml_m, 'ml_m', regressor=False, classifier=True)
-        self._learner = {'ml_m_z': clone(ml_m),
-                         'ml_g_du_z0': clone(ml_g), 'ml_g_du_z1': clone(ml_g),
-                         'ml_m_d_z0': clone(ml_m), 'ml_m_d_z1': clone(ml_m)}
-        self._predict_method = {'ml_m_z': 'predict_proba',
-                                'ml_g_du_z0': 'predict_proba', 'ml_g_du_z1': 'predict_proba',
-                                'ml_m_d_z0': 'predict_proba', 'ml_m_d_z1': 'predict_proba'}
+        _ = self._check_learner(ml_g, "ml_g", regressor=False, classifier=True)
+        _ = self._check_learner(ml_m, "ml_m", regressor=False, classifier=True)
+        self._learner = {
+            "ml_m_z": clone(ml_m),
+            "ml_g_du_z0": clone(ml_g),
+            "ml_g_du_z1": clone(ml_g),
+            "ml_m_d_z0": clone(ml_m),
+            "ml_m_d_z1": clone(ml_m),
+        }
+        self._predict_method = {
+            "ml_m_z": "predict_proba",
+            "ml_g_du_z0": "predict_proba",
+            "ml_g_du_z1": "predict_proba",
+            "ml_m_d_z0": "predict_proba",
+            "ml_m_d_z1": "predict_proba",
+        }
 
         self._initialize_ml_nuisance_params()
 
         if draw_sample_splitting:
             strata = self._dml_data.d.reshape(-1, 1) + 2 * self._dml_data.z.reshape(-1, 1)
-            obj_dml_resampling = DoubleMLResampling(n_folds=self.n_folds,
-                                                    n_rep=self.n_rep,
-                                                    n_obs=self._dml_data.n_obs,
-                                                    apply_cross_fitting=self.apply_cross_fitting,
-                                                    stratify=strata)
+            obj_dml_resampling = DoubleMLResampling(
+                n_folds=self.n_folds,
+                n_rep=self.n_rep,
+                n_obs=self._dml_data.n_obs,
+                apply_cross_fitting=self.apply_cross_fitting,
+                stratify=strata,
+            )
             self._smpls = obj_dml_resampling.split_samples()
+        self._external_predictions_implemented = True
 
     @property
     def quantile(self):
@@ -219,33 +233,33 @@ def trimming_threshold(self):
 
     @property
     def _score_element_names(self):
-        return ['ind_d', 'm_z', 'g_du_z0', 'g_du_z1', 'y', 'z', 'comp_prob']
+        return ["ind_d", "m_z", "g_du_z0", "g_du_z1", "y", "z", "comp_prob"]
 
     def _compute_ipw_score(self, theta, d, y, prop, z, comp_prob):
         sign = 2 * self.treatment - 1.0
         weights = sign * (z / prop - (1 - z) / (1 - prop)) / comp_prob
         u = (d == self._treatment) * (y <= theta)
-        v = -1. * self.quantile
+        v = -1.0 * self.quantile
         score = weights * u + v
         return score
 
     def _compute_score(self, psi_elements, coef, inds=None):
         sign = 2 * self.treatment - 1.0
-        ind_d = psi_elements['ind_d']
-        m_z = psi_elements['m_z']
-        g_du_z0 = psi_elements['g_du_z0']
-        g_du_z1 = psi_elements['g_du_z1']
-        y = psi_elements['y']
-        z = psi_elements['z']
-        comp_prob = psi_elements['comp_prob']
+        ind_d = psi_elements["ind_d"]
+        m_z = psi_elements["m_z"]
+        g_du_z0 = psi_elements["g_du_z0"]
+        g_du_z1 = psi_elements["g_du_z1"]
+        y = psi_elements["y"]
+        z = psi_elements["z"]
+        comp_prob = psi_elements["comp_prob"]
 
         if inds is not None:
-            ind_d = psi_elements['ind_d'][inds]
-            m_z = psi_elements['m_z']
-            g_du_z0 = psi_elements['g_du_z0'][inds]
-            g_du_z1 = psi_elements['g_du_z1'][inds]
-            y = psi_elements['y'][inds]
-            z = psi_elements['z'][inds]
+            ind_d = psi_elements["ind_d"][inds]
+            m_z = psi_elements["m_z"]
+            g_du_z0 = psi_elements["g_du_z0"][inds]
+            g_du_z1 = psi_elements["g_du_z1"][inds]
+            y = psi_elements["y"][inds]
+            z = psi_elements["z"][inds]
 
         score1 = g_du_z1 - g_du_z0
         score2 = (z / m_z) * (ind_d * (y <= coef) - g_du_z1)
@@ -255,17 +269,17 @@ def _compute_score(self, psi_elements, coef, inds=None):
 
     def _compute_score_deriv(self, psi_elements, coef, inds=None):
         sign = 2 * self.treatment - 1.0
-        ind_d = psi_elements['ind_d']
-        y = psi_elements['y']
-        m_z = psi_elements['m_z']
-        z = psi_elements['z']
-        comp_prob = psi_elements['comp_prob']
+        ind_d = psi_elements["ind_d"]
+        y = psi_elements["y"]
+        m_z = psi_elements["m_z"]
+        z = psi_elements["z"]
+        comp_prob = psi_elements["comp_prob"]
 
         if inds is not None:
-            ind_d = psi_elements['ind_d'][inds]
-            y = psi_elements['y'][inds]
-            m_z = psi_elements['m_z'][inds]
-            z = psi_elements['z'][inds]
+            ind_d = psi_elements["ind_d"][inds]
+            y = psi_elements["y"][inds]
+            m_z = psi_elements["m_z"][inds]
+            z = psi_elements["z"][inds]
 
         score_weights = sign * ((z / m_z) - (1 - z) / (1 - m_z)) * ind_d / comp_prob
         u = (y - coef).reshape(-1, 1)
@@ -274,179 +288,242 @@ def _compute_score_deriv(self, psi_elements, coef, inds=None):
         return deriv
 
     def _initialize_ml_nuisance_params(self):
-        self._params = {learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols}
-                        for learner in ['ml_m_z', 'ml_g_du_z0', 'ml_g_du_z1',
-                                        'ml_m_d_z0', 'ml_m_d_z1']}
-
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
-        x, y = check_X_y(self._dml_data.x, self._dml_data.y,
-                         force_all_finite=False)
-        x, d = check_X_y(x, self._dml_data.d,
-                         force_all_finite=False)
-        x, z = check_X_y(x, np.ravel(self._dml_data.z),
-                         force_all_finite=False)
+        self._params = {
+            learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols}
+            for learner in ["ml_m_z", "ml_g_du_z0", "ml_g_du_z1", "ml_m_d_z0", "ml_m_d_z1"]
+        }
+
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
+        x, y = check_X_y(self._dml_data.x, self._dml_data.y, force_all_finite=False)
+        x, d = check_X_y(x, self._dml_data.d, force_all_finite=False)
+        x, z = check_X_y(x, np.ravel(self._dml_data.z), force_all_finite=False)
+
+        m_z = external_predictions["ml_m_z"] is not None
+        m_d_d0 = external_predictions["ml_m_d_z0"] is not None
+        m_d_d1 = external_predictions["ml_m_d_z1"] is not None
+        g_du_z0 = external_predictions["ml_g_du_z0"] is not None
+        g_du_z1 = external_predictions["ml_g_du_z1"] is not None
+        ext_preds = [m_z, m_d_d0, m_d_d1, g_du_z0, g_du_z1]
 
         # create strata for splitting
         strata = self._dml_data.d.reshape(-1, 1) + 2 * self._dml_data.z.reshape(-1, 1)
 
         # initialize nuisance predictions, targets and models
-        m_z_hat = {'models': None,
-                   'targets': np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
-                   'preds': np.full(shape=self._dml_data.n_obs, fill_value=np.nan)
-                   }
-        m_d_z0_hat = copy.deepcopy(m_z_hat)
-        m_d_z1_hat = copy.deepcopy(m_z_hat)
-        g_du_z0_hat = copy.deepcopy(m_z_hat)
-        g_du_z1_hat = copy.deepcopy(m_z_hat)
-
-        # initialize models
-        fitted_models = {}
-        for learner in self.params_names:
-            # set nuisance model parameters
-            est_params = self._get_params(learner)
-            if est_params is not None:
-                fitted_models[learner] = [clone(self._learner[learner]).set_params(**est_params[i_fold])
-                                          for i_fold in range(self.n_folds)]
-            else:
-                fitted_models[learner] = [clone(self._learner[learner]) for i_fold in range(self.n_folds)]
+        if not all(ext_preds):
+            m_z_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+            m_d_z0_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+            m_d_z1_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+            g_du_z0_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+            g_du_z1_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+
+            # initialize models
+            fitted_models = {}
+            for learner in self.params_names:
+                # set nuisance model parameters
+                est_params = self._get_params(learner)
+                if est_params is not None:
+                    fitted_models[learner] = [
+                        clone(self._learner[learner]).set_params(**est_params[i_fold]) for i_fold in range(self.n_folds)
+                    ]
+                else:
+                    fitted_models[learner] = [clone(self._learner[learner]) for i_fold in range(self.n_folds)]
+            ipw_vec = np.full(shape=self.n_folds, fill_value=np.nan)
+        elif any(ext_preds) and not any(ext_preds):
+            raise ValueError("External predictions for all estimations or for none are required.")
+        else:
+            m_z_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_m_z"],
+            }
+            m_d_z0_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_m_d_z0"],
+            }
+            m_d_z1_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_m_d_z1"],
+            }
+            g_du_z0_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_g_du_z0"],
+            }
+            g_du_z1_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_g_du_z1"],
+            }
 
-        ipw_vec = np.full(shape=self.n_folds, fill_value=np.nan)
         # calculate nuisance functions over different folds
-        for i_fold in range(self.n_folds):
-            train_inds = smpls[i_fold][0]
-            test_inds = smpls[i_fold][1]
-
-            # start nested crossfitting
-            train_inds_1, train_inds_2 = train_test_split(train_inds, test_size=0.5,
-                                                          random_state=42, stratify=strata[train_inds])
-            smpls_prelim = [(train, test) for train, test in
-                            StratifiedKFold(n_splits=self.n_folds).split(X=train_inds_1, y=strata[train_inds_1])]
-
-            d_train_1 = d[train_inds_1]
-            y_train_1 = y[train_inds_1]
-            x_train_1 = x[train_inds_1, :]
-            z_train_1 = z[train_inds_1]
-
-            # preliminary propensity for z
-            ml_m_z_prelim = clone(fitted_models['ml_m_z'][i_fold])
-            m_z_hat_prelim = _dml_cv_predict(ml_m_z_prelim, x_train_1, z_train_1,
-                                             method='predict_proba', smpls=smpls_prelim)['preds']
-
-            m_z_hat_prelim = _trimm(m_z_hat_prelim, self.trimming_rule, self.trimming_threshold)
-            if self._normalize_ipw:
-                m_z_hat_prelim = _normalize_ipw(m_z_hat_prelim, z_train_1)
-
-            # propensity for d == 1 cond. on z == 0 (training set 1)
-            z0_train_1 = z_train_1 == 0
-            x_z0_train_1 = x_train_1[z0_train_1, :]
-            d_z0_train_1 = d_train_1[z0_train_1]
-            ml_m_d_z0_prelim = clone(fitted_models['ml_m_d_z0'][i_fold])
-            ml_m_d_z0_prelim.fit(x_z0_train_1, d_z0_train_1)
-            m_d_z0_hat_prelim = _predict_zero_one_propensity(ml_m_d_z0_prelim, x_train_1)
-
-            # propensity for d == 1 cond. on z == 1 (training set 1)
-            z1_train_1 = z_train_1 == 1
-            x_z1_train_1 = x_train_1[z1_train_1, :]
-            d_z1_train_1 = d_train_1[z1_train_1]
-            ml_m_d_z1_prelim = clone(fitted_models['ml_m_d_z1'][i_fold])
-            ml_m_d_z1_prelim.fit(x_z1_train_1, d_z1_train_1)
-            m_d_z1_hat_prelim = _predict_zero_one_propensity(ml_m_d_z1_prelim, x_train_1)
-
-            # preliminary estimate of theta_2_aux
-            comp_prob_prelim = np.mean(m_d_z1_hat_prelim - m_d_z0_hat_prelim
-                                       + z_train_1 / m_z_hat_prelim * (d_train_1 - m_d_z1_hat_prelim)
-                                       - (1 - z_train_1) / (1 - m_z_hat_prelim) * (d_train_1 - m_d_z0_hat_prelim))
-
-            # preliminary ipw estimate
-            def ipw_score(theta):
-                res = np.mean(self._compute_ipw_score(theta, d_train_1, y_train_1, m_z_hat_prelim,
-                                                      z_train_1, comp_prob_prelim))
-                return res
-
-            _, bracket_guess = _get_bracket_guess(ipw_score, self._coef_start_val, self._coef_bounds)
-            ipw_est = _solve_ipw_score(ipw_score=ipw_score, bracket_guess=bracket_guess)
-            ipw_vec[i_fold] = ipw_est
-
-            # use the preliminary estimates to fit the nuisance parameters on train_2
-            d_train_2 = d[train_inds_2]
-            y_train_2 = y[train_inds_2]
-            x_train_2 = x[train_inds_2, :]
-            z_train_2 = z[train_inds_2]
-
-            # define test observations
-            d_test = d[test_inds]
-            y_test = y[test_inds]
-            x_test = x[test_inds, :]
-            z_test = z[test_inds]
-
-            # propensity for (D == treatment)*Ind(Y <= ipq_est) cond. on z == 0
-            z0_train_2 = z_train_2 == 0
-            x_z0_train_2 = x_train_2[z0_train_2, :]
-            du_z0_train_2 = (d_train_2[z0_train_2] == self._treatment) * (y_train_2[z0_train_2] <= ipw_est)
-            fitted_models['ml_g_du_z0'][i_fold].fit(x_z0_train_2, du_z0_train_2)
-            g_du_z0_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_g_du_z0'][i_fold], x_test)
-
-            # propensity for (D == treatment)*Ind(Y <= ipq_est) cond. on z == 1
-            z1_train_2 = z_train_2 == 1
-            x_z1_train_2 = x_train_2[z1_train_2, :]
-            du_z1_train_2 = (d_train_2[z1_train_2] == self._treatment) * (y_train_2[z1_train_2] <= ipw_est)
-            fitted_models['ml_g_du_z1'][i_fold].fit(x_z1_train_2, du_z1_train_2)
-            g_du_z1_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_g_du_z1'][i_fold], x_test)
-
-            # the predictions of both should only be evaluated conditional on z == 0 or z == 1
-            test_inds_z0 = test_inds[z_test == 0]
-            test_inds_z1 = test_inds[z_test == 1]
-            g_du_z0_hat['targets'][test_inds_z0] = (1.0 * (d_test[z_test == 0] == self._treatment) *
-                                                    (y_test[z_test == 0] <= ipw_est))
-            g_du_z1_hat['targets'][test_inds_z1] = (1.0 * (d_test[z_test == 1] == self._treatment) *
-                                                    (y_test[z_test == 1] <= ipw_est))
-
-            # refit nuisance elements for the local potential quantile
-            z_train = z[train_inds]
-            x_train = x[train_inds]
-            d_train = d[train_inds]
-
-            # refit propensity for z (whole training set)
-            fitted_models['ml_m_z'][i_fold].fit(x_train, z_train)
-            m_z_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_m_z'][i_fold], x_test)
-
-            # refit propensity for d == 1 cond. on z == 0 (whole training set)
-            z0_train = z_train == 0
-            x_z0_train = x_train[z0_train, :]
-            d_z0_train = d_train[z0_train]
-            fitted_models['ml_m_d_z0'][i_fold].fit(x_z0_train, d_z0_train)
-            m_d_z0_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_m_d_z0'][i_fold], x_test)
-
-            # propensity for d == 1 cond. on z == 1 (whole training set)
-            x_z1_train = x_train[z_train == 1, :]
-            d_z1_train = d_train[z_train == 1]
-            fitted_models['ml_m_d_z1'][i_fold].fit(x_z1_train, d_z1_train)
-            m_d_z1_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_m_d_z1'][i_fold], x_test)
+        if not all(ext_preds):
+            for i_fold in range(self.n_folds):
+                train_inds = smpls[i_fold][0]
+                test_inds = smpls[i_fold][1]
+
+                # start nested crossfitting
+                train_inds_1, train_inds_2 = train_test_split(
+                    train_inds, test_size=0.5, random_state=42, stratify=strata[train_inds]
+                )
+                smpls_prelim = [
+                    (train, test)
+                    for train, test in StratifiedKFold(n_splits=self.n_folds).split(X=train_inds_1, y=strata[train_inds_1])
+                ]
+
+                d_train_1 = d[train_inds_1]
+                y_train_1 = y[train_inds_1]
+                x_train_1 = x[train_inds_1, :]
+                z_train_1 = z[train_inds_1]
+
+                # preliminary propensity for z
+                ml_m_z_prelim = clone(fitted_models["ml_m_z"][i_fold])
+                m_z_hat_prelim = _dml_cv_predict(ml_m_z_prelim, x_train_1, z_train_1,
+                                                 method="predict_proba", smpls=smpls_prelim)[
+                    "preds"
+                ]
+
+                m_z_hat_prelim = _trimm(m_z_hat_prelim, self.trimming_rule, self.trimming_threshold)
+                if self._normalize_ipw:
+                    m_z_hat_prelim = _normalize_ipw(m_z_hat_prelim, z_train_1)
+
+                # propensity for d == 1 cond. on z == 0 (training set 1)
+                z0_train_1 = z_train_1 == 0
+                x_z0_train_1 = x_train_1[z0_train_1, :]
+                d_z0_train_1 = d_train_1[z0_train_1]
+                ml_m_d_z0_prelim = clone(fitted_models["ml_m_d_z0"][i_fold])
+                ml_m_d_z0_prelim.fit(x_z0_train_1, d_z0_train_1)
+                m_d_z0_hat_prelim = _predict_zero_one_propensity(ml_m_d_z0_prelim, x_train_1)
+
+                # propensity for d == 1 cond. on z == 1 (training set 1)
+                z1_train_1 = z_train_1 == 1
+                x_z1_train_1 = x_train_1[z1_train_1, :]
+                d_z1_train_1 = d_train_1[z1_train_1]
+                ml_m_d_z1_prelim = clone(fitted_models["ml_m_d_z1"][i_fold])
+                ml_m_d_z1_prelim.fit(x_z1_train_1, d_z1_train_1)
+                m_d_z1_hat_prelim = _predict_zero_one_propensity(ml_m_d_z1_prelim, x_train_1)
+
+                # preliminary estimate of theta_2_aux
+                comp_prob_prelim = np.mean(
+                    m_d_z1_hat_prelim
+                    - m_d_z0_hat_prelim
+                    + z_train_1 / m_z_hat_prelim * (d_train_1 - m_d_z1_hat_prelim)
+                    - (1 - z_train_1) / (1 - m_z_hat_prelim) * (d_train_1 - m_d_z0_hat_prelim)
+                )
+
+                # preliminary ipw estimate
+                def ipw_score(theta):
+                    res = np.mean(
+                        self._compute_ipw_score(theta, d_train_1, y_train_1, m_z_hat_prelim, z_train_1, comp_prob_prelim)
+                    )
+                    return res
+
+                _, bracket_guess = _get_bracket_guess(ipw_score, self._coef_start_val, self._coef_bounds)
+                ipw_est = _solve_ipw_score(ipw_score=ipw_score, bracket_guess=bracket_guess)
+                ipw_vec[i_fold] = ipw_est
+
+                # use the preliminary estimates to fit the nuisance parameters on train_2
+                d_train_2 = d[train_inds_2]
+                y_train_2 = y[train_inds_2]
+                x_train_2 = x[train_inds_2, :]
+                z_train_2 = z[train_inds_2]
+
+                # define test observations
+                d_test = d[test_inds]
+                y_test = y[test_inds]
+                x_test = x[test_inds, :]
+                z_test = z[test_inds]
+
+                # propensity for (D == treatment)*Ind(Y <= ipq_est) cond. on z == 0
+                z0_train_2 = z_train_2 == 0
+                x_z0_train_2 = x_train_2[z0_train_2, :]
+                du_z0_train_2 = (d_train_2[z0_train_2] == self._treatment) * (y_train_2[z0_train_2] <= ipw_est)
+                fitted_models["ml_g_du_z0"][i_fold].fit(x_z0_train_2, du_z0_train_2)
+                g_du_z0_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_g_du_z0"][i_fold], x_test)
+
+                # propensity for (D == treatment)*Ind(Y <= ipq_est) cond. on z == 1
+                z1_train_2 = z_train_2 == 1
+                x_z1_train_2 = x_train_2[z1_train_2, :]
+                du_z1_train_2 = (d_train_2[z1_train_2] == self._treatment) * (y_train_2[z1_train_2] <= ipw_est)
+                fitted_models["ml_g_du_z1"][i_fold].fit(x_z1_train_2, du_z1_train_2)
+                g_du_z1_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_g_du_z1"][i_fold], x_test)
+
+                # the predictions of both should only be evaluated conditional on z == 0 or z == 1
+                test_inds_z0 = test_inds[z_test == 0]
+                test_inds_z1 = test_inds[z_test == 1]
+                g_du_z0_hat["targets"][test_inds_z0] = (
+                    1.0 * (d_test[z_test == 0] == self._treatment) * (y_test[z_test == 0] <= ipw_est)
+                )
+                g_du_z1_hat["targets"][test_inds_z1] = (
+                    1.0 * (d_test[z_test == 1] == self._treatment) * (y_test[z_test == 1] <= ipw_est)
+                )
+
+                # refit nuisance elements for the local potential quantile
+                z_train = z[train_inds]
+                x_train = x[train_inds]
+                d_train = d[train_inds]
+
+                # refit propensity for z (whole training set)
+                fitted_models["ml_m_z"][i_fold].fit(x_train, z_train)
+                m_z_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_m_z"][i_fold], x_test)
+
+                # refit propensity for d == 1 cond. on z == 0 (whole training set)
+                z0_train = z_train == 0
+                x_z0_train = x_train[z0_train, :]
+                d_z0_train = d_train[z0_train]
+                fitted_models["ml_m_d_z0"][i_fold].fit(x_z0_train, d_z0_train)
+                m_d_z0_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_m_d_z0"][i_fold], x_test)
+
+                # propensity for d == 1 cond. on z == 1 (whole training set)
+                x_z1_train = x_train[z_train == 1, :]
+                d_z1_train = d_train[z_train == 1]
+                fitted_models["ml_m_d_z1"][i_fold].fit(x_z1_train, d_z1_train)
+                m_d_z1_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_m_d_z1"][i_fold], x_test)
 
         # save targets and models
-        m_z_hat['targets'] = z
+        m_z_hat["targets"] = z
+
         # set targets to relevant subsample
-        g_du_z0_hat['targets'][z == 1] = np.nan
-        g_du_z1_hat['targets'][z == 0] = np.nan
+        g_du_z0_hat["targets"] = _cond_targets(g_du_z0_hat["targets"], cond_sample=(z == 0))
+        g_du_z1_hat["targets"] = _cond_targets(g_du_z1_hat["targets"], cond_sample=(z == 1))
 
         # the predictions of both should only be evaluated conditional on z == 0 or z == 1
-        m_d_z0_hat['targets'][z == 0] = d[z == 0]
-        m_d_z0_hat['targets'][z == 1] = np.nan
-        m_d_z1_hat['targets'][z == 1] = d[z == 1]
-        m_d_z1_hat['targets'][z == 0] = np.nan
+        m_d_z0_hat["targets"] = _cond_targets(d, cond_sample=(z == 0))
+        m_d_z0_hat["targets"] = _cond_targets(d, cond_sample=(z == 1))
 
         if return_models:
-            m_z_hat['models'] = fitted_models['ml_m_z']
-            m_d_z0_hat['models'] = fitted_models['ml_m_d_z0']
-            m_d_z1_hat['models'] = fitted_models['ml_m_d_z1']
-            g_du_z0_hat['models'] = fitted_models['ml_g_du_z0']
-            g_du_z1_hat['models'] = fitted_models['ml_g_du_z1']
+            m_z_hat["models"] = fitted_models["ml_m_z"]
+            m_d_z0_hat["models"] = fitted_models["ml_m_d_z0"]
+            m_d_z1_hat["models"] = fitted_models["ml_m_d_z1"]
+            g_du_z0_hat["models"] = fitted_models["ml_g_du_z0"]
+            g_du_z1_hat["models"] = fitted_models["ml_g_du_z1"]
 
         # clip propensities
-        m_z_hat_adj = _trimm(m_z_hat['preds'], self.trimming_rule, self.trimming_threshold)
+        m_z_hat_adj = _trimm(m_z_hat["preds"], self.trimming_rule, self.trimming_threshold)
 
         if self._normalize_ipw:
-            if self.dml_procedure == 'dml1':
+            if self.dml_procedure == "dml1":
                 for _, test_index in smpls:
                     m_z_hat_adj[test_index] = _normalize_ipw(m_z_hat_adj[test_index], z[test_index])
             else:
@@ -454,49 +531,60 @@ def ipw_score(theta):
 
         # this could be adjusted to be compatible with dml1
         # estimate final nuisance parameter
-        comp_prob_hat = np.mean(m_d_z1_hat['preds'] - m_d_z0_hat['preds']
-                                + z / m_z_hat_adj * (d - m_d_z1_hat['preds'])
-                                - (1 - z) / (1 - m_z_hat_adj) * (d - m_d_z0_hat['preds']))
-
-        # readjust start value for minimization
-        self._coef_start_val = np.mean(ipw_vec)
-
-        psi_elements = {'ind_d': d == self._treatment, 'm_z': m_z_hat_adj,
-                        'g_du_z0': g_du_z0_hat['preds'], 'g_du_z1': g_du_z1_hat['preds'],
-                        'y': y, 'z': z, 'comp_prob': comp_prob_hat}
-        preds = {'predictions': {'ml_m_z':  m_z_hat['preds'],
-                                 'ml_m_d_z0': m_d_z0_hat['preds'],
-                                 'ml_m_d_z1': m_d_z1_hat['preds'],
-                                 'ml_g_du_z0': g_du_z0_hat['preds'],
-                                 'ml_g_du_z1': g_du_z1_hat['preds']},
-                 'targets': {'ml_m_z':  m_z_hat['targets'],
-                             'ml_m_d_z0': m_d_z0_hat['targets'],
-                             'ml_m_d_z1': m_d_z1_hat['targets'],
-                             'ml_g_du_z0': g_du_z0_hat['targets'],
-                             'ml_g_du_z1': g_du_z1_hat['targets']},
-                 'models': {'ml_m_z':  m_z_hat['models'],
-                            'ml_m_d_z0': m_d_z0_hat['models'],
-                            'ml_m_d_z1': m_d_z1_hat['models'],
-                            'ml_g_du_z0': g_du_z0_hat['models'],
-                            'ml_g_du_z1': g_du_z1_hat['models']}
-                 }
+        comp_prob_hat = np.mean(
+            m_d_z1_hat["preds"]
+            - m_d_z0_hat["preds"]
+            + z / m_z_hat_adj * (d - m_d_z1_hat["preds"])
+            - (1 - z) / (1 - m_z_hat_adj) * (d - m_d_z0_hat["preds"])
+        )
+
+        if not all(ext_preds):
+            # readjust start value for minimization
+            self._coef_start_val = np.mean(ipw_vec)
+
+        psi_elements = {
+            "ind_d": d == self._treatment,
+            "m_z": m_z_hat_adj,
+            "g_du_z0": g_du_z0_hat["preds"],
+            "g_du_z1": g_du_z1_hat["preds"],
+            "y": y,
+            "z": z,
+            "comp_prob": comp_prob_hat,
+        }
+        preds = {
+            "predictions": {
+                "ml_m_z": m_z_hat["preds"],
+                "ml_m_d_z0": m_d_z0_hat["preds"],
+                "ml_m_d_z1": m_d_z1_hat["preds"],
+                "ml_g_du_z0": g_du_z0_hat["preds"],
+                "ml_g_du_z1": g_du_z1_hat["preds"],
+            },
+            "targets": {
+                "ml_m_z": m_z_hat["targets"],
+                "ml_m_d_z0": m_d_z0_hat["targets"],
+                "ml_m_d_z1": m_d_z1_hat["targets"],
+                "ml_g_du_z0": g_du_z0_hat["targets"],
+                "ml_g_du_z1": g_du_z1_hat["targets"],
+            },
+            "models": {
+                "ml_m_z": m_z_hat["models"],
+                "ml_m_d_z0": m_d_z0_hat["models"],
+                "ml_m_d_z1": m_d_z1_hat["models"],
+                "ml_g_du_z0": g_du_z0_hat["models"],
+                "ml_g_du_z1": g_du_z1_hat["models"],
+            },
+        }
         return psi_elements, preds
 
-    def _nuisance_tuning(self, smpls, param_grids, scoring_methods, n_folds_tune, n_jobs_cv,
-                         search_mode, n_iter_randomized_search):
-        x, y = check_X_y(self._dml_data.x, self._dml_data.y,
-                         force_all_finite=False)
-        x, d = check_X_y(x, self._dml_data.d,
-                         force_all_finite=False)
-        x, z = check_X_y(x, np.ravel(self._dml_data.z),
-                         force_all_finite=False)
+    def _nuisance_tuning(
+        self, smpls, param_grids, scoring_methods, n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search
+    ):
+        x, y = check_X_y(self._dml_data.x, self._dml_data.y, force_all_finite=False)
+        x, d = check_X_y(x, self._dml_data.d, force_all_finite=False)
+        x, z = check_X_y(x, np.ravel(self._dml_data.z), force_all_finite=False)
 
         if scoring_methods is None:
-            scoring_methods = {'ml_m_z': None,
-                               'ml_m_d_z0': None,
-                               'ml_m_d_z1': None,
-                               'ml_g_du_z0': None,
-                               'ml_g_du_z1': None}
+            scoring_methods = {"ml_m_z": None, "ml_m_d_z0": None, "ml_m_d_z1": None, "ml_g_du_z0": None, "ml_g_du_z1": None}
 
         train_inds = [train_index for (train_index, _) in smpls]
         train_inds_z0 = [np.intersect1d(np.where(z == 0)[0], train) for train, _ in smpls]
@@ -505,21 +593,66 @@ def _nuisance_tuning(self, smpls, param_grids, scoring_methods, n_folds_tune, n_
         approx_quant = np.quantile(y[d == self.treatment], self.quantile)
         du = (d == self.treatment) * (y <= approx_quant)
 
-        m_z_tune_res = _dml_tune(z, x, train_inds,
-                                 self._learner['ml_m_z'], param_grids['ml_m_z'], scoring_methods['ml_m_z'],
-                                 n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
-        m_d_z0_tune_res = _dml_tune(d, x, train_inds_z0,
-                                    self._learner['ml_m_d_z0'], param_grids['ml_m_d_z0'], scoring_methods['ml_m_d_z0'],
-                                    n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
-        m_d_z1_tune_res = _dml_tune(d, x, train_inds_z1,
-                                    self._learner['ml_m_d_z1'], param_grids['ml_m_d_z1'], scoring_methods['ml_m_d_z1'],
-                                    n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
-        g_du_z0_tune_res = _dml_tune(du, x, train_inds_z0,
-                                     self._learner['ml_g_du_z0'], param_grids['ml_g_du_z0'], scoring_methods['ml_g_du_z0'],
-                                     n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
-        g_du_z1_tune_res = _dml_tune(du, x, train_inds_z1,
-                                     self._learner['ml_g_du_z1'], param_grids['ml_g_du_z1'], scoring_methods['ml_g_du_z1'],
-                                     n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
+        m_z_tune_res = _dml_tune(
+            z,
+            x,
+            train_inds,
+            self._learner["ml_m_z"],
+            param_grids["ml_m_z"],
+            scoring_methods["ml_m_z"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
+        m_d_z0_tune_res = _dml_tune(
+            d,
+            x,
+            train_inds_z0,
+            self._learner["ml_m_d_z0"],
+            param_grids["ml_m_d_z0"],
+            scoring_methods["ml_m_d_z0"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
+        m_d_z1_tune_res = _dml_tune(
+            d,
+            x,
+            train_inds_z1,
+            self._learner["ml_m_d_z1"],
+            param_grids["ml_m_d_z1"],
+            scoring_methods["ml_m_d_z1"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
+        g_du_z0_tune_res = _dml_tune(
+            du,
+            x,
+            train_inds_z0,
+            self._learner["ml_g_du_z0"],
+            param_grids["ml_g_du_z0"],
+            scoring_methods["ml_g_du_z0"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
+        g_du_z1_tune_res = _dml_tune(
+            du,
+            x,
+            train_inds_z1,
+            self._learner["ml_g_du_z1"],
+            param_grids["ml_g_du_z1"],
+            scoring_methods["ml_g_du_z1"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
 
         m_z_best_params = [xx.best_params_ for xx in m_z_tune_res]
         m_d_z0_best_params = [xx.best_params_ for xx in m_d_z0_tune_res]
@@ -527,34 +660,40 @@ def _nuisance_tuning(self, smpls, param_grids, scoring_methods, n_folds_tune, n_
         g_du_z0_best_params = [xx.best_params_ for xx in g_du_z0_tune_res]
         g_du_z1_best_params = [xx.best_params_ for xx in g_du_z1_tune_res]
 
-        params = {'ml_m_z': m_z_best_params,
-                  'ml_m_d_z0': m_d_z0_best_params,
-                  'ml_m_d_z1': m_d_z1_best_params,
-                  'ml_g_du_z0': g_du_z0_best_params,
-                  'ml_g_du_z1': g_du_z1_best_params}
-        tune_res = {'ml_m_z': m_z_tune_res,
-                    'ml_m_d_z0': m_d_z0_tune_res,
-                    'ml_m_d_z1':  m_d_z1_tune_res,
-                    'ml_g_du_z0': g_du_z0_tune_res,
-                    'ml_g_du_z1': g_du_z1_tune_res}
-
-        res = {'params': params,
-               'tune_res': tune_res}
+        params = {
+            "ml_m_z": m_z_best_params,
+            "ml_m_d_z0": m_d_z0_best_params,
+            "ml_m_d_z1": m_d_z1_best_params,
+            "ml_g_du_z0": g_du_z0_best_params,
+            "ml_g_du_z1": g_du_z1_best_params,
+        }
+        tune_res = {
+            "ml_m_z": m_z_tune_res,
+            "ml_m_d_z0": m_d_z0_tune_res,
+            "ml_m_d_z1": m_d_z1_tune_res,
+            "ml_g_du_z0": g_du_z0_tune_res,
+            "ml_g_du_z1": g_du_z1_tune_res,
+        }
+
+        res = {"params": params, "tune_res": tune_res}
 
         return res
 
     def _check_data(self, obj_dml_data):
         if not isinstance(obj_dml_data, DoubleMLData):
-            raise TypeError('The data must be of DoubleMLData type. '
-                            f'{str(obj_dml_data)} of type {str(type(obj_dml_data))} was passed.')
+            raise TypeError(
+                "The data must be of DoubleMLData type. " f"{str(obj_dml_data)} of type {str(type(obj_dml_data))} was passed."
+            )
         _check_zero_one_treatment(self)
-        one_instr = (obj_dml_data.n_instr == 1)
-        err_msg = ('Incompatible data. '
-                   'To fit an LPQ model with DML '
-                   'exactly one binary variable with values 0 and 1 '
-                   'needs to be specified as instrumental variable.')
+        one_instr = obj_dml_data.n_instr == 1
+        err_msg = (
+            "Incompatible data. "
+            "To fit an LPQ model with DML "
+            "exactly one binary variable with values 0 and 1 "
+            "needs to be specified as instrumental variable."
+        )
         if one_instr:
-            binary_instr = (type_of_target(obj_dml_data.z) == 'binary')
+            binary_instr = type_of_target(obj_dml_data.z) == "binary"
             zero_one_instr = np.all((np.power(obj_dml_data.z, 2) - obj_dml_data.z) == 0)
             if not (one_instr & binary_instr & zero_one_instr):
                 raise ValueError(err_msg)
diff --git a/doubleml/double_ml_pliv.py b/doubleml/double_ml_pliv.py
index 4c70c495..6725e925 100644
--- a/doubleml/double_ml_pliv.py
+++ b/doubleml/double_ml_pliv.py
@@ -145,6 +145,7 @@ def __init__(self,
         if 'ml_g' in self._learner:
             self._predict_method['ml_g'] = 'predict'
         self._initialize_ml_nuisance_params()
+        self._external_predictions_implemented = True
 
     @classmethod
     def _partialX(cls,
@@ -283,9 +284,9 @@ def _check_data(self, obj_dml_data):
                              'use DoubleMLPLR instead of DoubleMLPLIV.')
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         if self.partialX & (not self.partialZ):
-            psi_elements, preds = self._nuisance_est_partial_x(smpls, n_jobs_cv, return_models)
+            psi_elements, preds = self._nuisance_est_partial_x(smpls, n_jobs_cv, external_predictions, return_models)
         elif (not self.partialX) & self.partialZ:
             psi_elements, preds = self._nuisance_est_partial_z(smpls, n_jobs_cv, return_models)
         else:
@@ -309,16 +310,21 @@ def _nuisance_tuning(self, smpls, param_grids, scoring_methods, n_folds_tune, n_
 
         return res
 
-    def _nuisance_est_partial_x(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est_partial_x(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
                          force_all_finite=False)
 
         # nuisance l
-        l_hat = _dml_cv_predict(self._learner['ml_l'], x, y, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_l'), method=self._predict_method['ml_l'],
-                                return_models=return_models)
+        if external_predictions['ml_l'] is not None:
+            l_hat = {'preds': external_predictions['ml_l'],
+                     'targets': None,
+                     'models': None}
+        else:
+            l_hat = _dml_cv_predict(self._learner['ml_l'], x, y, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_l'), method=self._predict_method['ml_l'],
+                                    return_models=return_models)
         _check_finite_predictions(l_hat['preds'], self._learner['ml_l'], 'ml_l', smpls)
 
         predictions = {'ml_l': l_hat['preds']}
@@ -329,36 +335,54 @@ def _nuisance_est_partial_x(self, smpls, n_jobs_cv, return_models=False):
             # one instrument: just identified
             x, z = check_X_y(x, np.ravel(self._dml_data.z),
                              force_all_finite=False)
-            m_hat = _dml_cv_predict(self._learner['ml_m'], x, z, smpls=smpls, n_jobs=n_jobs_cv,
-                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                    return_models=return_models)
+            if external_predictions['ml_m'] is not None:
+                m_hat = {'preds': external_predictions['ml_m'],
+                         'targets': None,
+                         'models': None}
+            else:
+                m_hat = _dml_cv_predict(self._learner['ml_m'], x, z, smpls=smpls, n_jobs=n_jobs_cv,
+                                        est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                        return_models=return_models)
             predictions['ml_m'] = m_hat['preds']
             targets['ml_m'] = m_hat['targets']
             models['ml_m'] = m_hat['models']
         else:
             # several instruments: 2SLS
             m_hat = {'preds': np.full((self._dml_data.n_obs, self._dml_data.n_instr), np.nan),
+                     'targets': [None] * self._dml_data.n_instr,
                      'models': [None] * self._dml_data.n_instr}
-            z = self._dml_data.z
             for i_instr in range(self._dml_data.n_instr):
+                z = self._dml_data.z
                 x, this_z = check_X_y(x, z[:, i_instr],
                                       force_all_finite=False)
-                res_cv_predict = _dml_cv_predict(self._learner['ml_m'], x, this_z, smpls=smpls, n_jobs=n_jobs_cv,
-                                                 est_params=self._get_params('ml_m_' + self._dml_data.z_cols[i_instr]),
-                                                 method=self._predict_method['ml_m'], return_models=return_models)
+                if external_predictions['ml_m_' + self._dml_data.z_cols[i_instr]] is not None:
+                    m_hat['preds'][:, i_instr] = external_predictions['ml_m_' + self._dml_data.z_cols[i_instr]]
+                    predictions['ml_m_' + self._dml_data.z_cols[i_instr]] = external_predictions[
+                        'ml_m_' + self._dml_data.z_cols[i_instr]]
+                    targets['ml_m_' + self._dml_data.z_cols[i_instr]] = None
+                    models['ml_m_' + self._dml_data.z_cols[i_instr]] = None
+                else:
+                    res_cv_predict = _dml_cv_predict(self._learner['ml_m'], x, this_z, smpls=smpls, n_jobs=n_jobs_cv,
+                                                     est_params=self._get_params('ml_m_' + self._dml_data.z_cols[i_instr]),
+                                                     method=self._predict_method['ml_m'], return_models=return_models)
 
-                m_hat['preds'][:, i_instr] = res_cv_predict['preds']
+                    m_hat['preds'][:, i_instr] = res_cv_predict['preds']
 
-                predictions['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['preds']
-                targets['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['targets']
-                models['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['models']
+                    predictions['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['preds']
+                    targets['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['targets']
+                    models['ml_m_' + self._dml_data.z_cols[i_instr]] = res_cv_predict['models']
 
         _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
 
         # nuisance r
-        r_hat = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_r'), method=self._predict_method['ml_r'],
-                                return_models=return_models)
+        if external_predictions['ml_r'] is not None:
+            r_hat = {'preds': external_predictions['ml_r'],
+                     'targets': None,
+                     'models': None}
+        else:
+            r_hat = _dml_cv_predict(self._learner['ml_r'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_r'), method=self._predict_method['ml_r'],
+                                    return_models=return_models)
         _check_finite_predictions(r_hat['preds'], self._learner['ml_r'], 'ml_r', smpls)
         predictions['ml_r'] = r_hat['preds']
         targets['ml_r'] = r_hat['targets']
@@ -368,13 +392,18 @@ def _nuisance_est_partial_x(self, smpls, n_jobs_cv, return_models=False):
         if (self._dml_data.n_instr == 1) & ('ml_g' in self._learner):
             # an estimate of g is obtained for the IV-type score and callable scores
             # get an initial estimate for theta using the partialling out score
-            psi_a = -np.multiply(d - r_hat['preds'], z - m_hat['preds'])
-            psi_b = np.multiply(z - m_hat['preds'], y - l_hat['preds'])
-            theta_initial = -np.nanmean(psi_b) / np.nanmean(psi_a)
-            # nuisance g
-            g_hat = _dml_cv_predict(self._learner['ml_g'], x, y - theta_initial * d, smpls=smpls, n_jobs=n_jobs_cv,
-                                    est_params=self._get_params('ml_g'), method=self._predict_method['ml_g'],
-                                    return_models=return_models)
+            if external_predictions['ml_g'] is not None:
+                g_hat = {'preds': external_predictions['ml_g'],
+                         'targets': None,
+                         'models': None}
+            else:
+                psi_a = -np.multiply(d - r_hat['preds'], z - m_hat['preds'])
+                psi_b = np.multiply(z - m_hat['preds'], y - l_hat['preds'])
+                theta_initial = -np.nanmean(psi_b) / np.nanmean(psi_a)
+                # nuisance g
+                g_hat = _dml_cv_predict(self._learner['ml_g'], x, y - theta_initial * d, smpls=smpls, n_jobs=n_jobs_cv,
+                                        est_params=self._get_params('ml_g'), method=self._predict_method['ml_g'],
+                                        return_models=return_models)
             _check_finite_predictions(g_hat['preds'], self._learner['ml_g'], 'ml_g', smpls)
 
         predictions['ml_g'] = g_hat['preds']
diff --git a/doubleml/double_ml_plr.py b/doubleml/double_ml_plr.py
index d75cbcf8..41e7df2f 100644
--- a/doubleml/double_ml_plr.py
+++ b/doubleml/double_ml_plr.py
@@ -150,6 +150,7 @@ def __init__(self,
 
         self._initialize_ml_nuisance_params()
         self._sensitivity_implemented = True
+        self._external_predictions_implemented = True
 
     def _initialize_ml_nuisance_params(self):
         self._params = {learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols}
@@ -166,23 +167,39 @@ def _check_data(self, obj_dml_data):
                              'To fit a partially linear IV regression model use DoubleMLPLIV instead of DoubleMLPLR.')
         return
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
                          force_all_finite=False)
+        m_external = external_predictions['ml_m'] is not None
+        l_external = external_predictions['ml_l'] is not None
+        if 'ml_g' in self._learner:
+            g_external = external_predictions['ml_g'] is not None
+        else:
+            g_external = False
 
         # nuisance l
-        l_hat = _dml_cv_predict(self._learner['ml_l'], x, y, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_l'), method=self._predict_method['ml_l'],
-                                return_models=return_models)
-        _check_finite_predictions(l_hat['preds'], self._learner['ml_l'], 'ml_l', smpls)
+        if l_external:
+            l_hat = {'preds': external_predictions['ml_l'],
+                     'targets': None,
+                     'models': None}
+        else:
+            l_hat = _dml_cv_predict(self._learner['ml_l'], x, y, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_l'), method=self._predict_method['ml_l'],
+                                    return_models=return_models)
+            _check_finite_predictions(l_hat['preds'], self._learner['ml_l'], 'ml_l', smpls)
 
         # nuisance m
-        m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
-                                est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
-                                return_models=return_models)
-        _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
+        if m_external:
+            m_hat = {'preds': external_predictions['ml_m'],
+                     'targets': None,
+                     'models': None}
+        else:
+            m_hat = _dml_cv_predict(self._learner['ml_m'], x, d, smpls=smpls, n_jobs=n_jobs_cv,
+                                    est_params=self._get_params('ml_m'), method=self._predict_method['ml_m'],
+                                    return_models=return_models)
+            _check_finite_predictions(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls)
         if self._check_learner(self._learner['ml_m'], 'ml_m', regressor=True, classifier=True):
             _check_is_propensity(m_hat['preds'], self._learner['ml_m'], 'ml_m', smpls, eps=1e-12)
 
@@ -203,10 +220,15 @@ def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
             psi_b = np.multiply(d - m_hat['preds'], y - l_hat['preds'])
             theta_initial = -np.nanmean(psi_b) / np.nanmean(psi_a)
             # nuisance g
-            g_hat = _dml_cv_predict(self._learner['ml_g'], x, y - theta_initial*d, smpls=smpls, n_jobs=n_jobs_cv,
-                                    est_params=self._get_params('ml_g'), method=self._predict_method['ml_g'],
-                                    return_models=return_models)
-            _check_finite_predictions(g_hat['preds'], self._learner['ml_g'], 'ml_g', smpls)
+            if g_external:
+                g_hat = {'preds': external_predictions['ml_g'],
+                         'targets': None,
+                         'models': None}
+            else:
+                g_hat = _dml_cv_predict(self._learner['ml_g'], x, y - theta_initial*d, smpls=smpls, n_jobs=n_jobs_cv,
+                                        est_params=self._get_params('ml_g'), method=self._predict_method['ml_g'],
+                                        return_models=return_models)
+                _check_finite_predictions(g_hat['preds'], self._learner['ml_g'], 'ml_g', smpls)
 
         psi_a, psi_b = self._score_elements(y, d, l_hat['preds'], m_hat['preds'], g_hat['preds'], smpls)
         psi_elements = {'psi_a': psi_a,
diff --git a/doubleml/double_ml_pq.py b/doubleml/double_ml_pq.py
index d9a36d41..dee3dc12 100644
--- a/doubleml/double_ml_pq.py
+++ b/doubleml/double_ml_pq.py
@@ -1,5 +1,4 @@
 import numpy as np
-import copy
 from sklearn.base import clone
 from sklearn.utils import check_X_y
 from sklearn.model_selection import StratifiedKFold, train_test_split
@@ -8,11 +7,26 @@
 from .double_ml_score_mixins import NonLinearScoreMixin
 from .double_ml_data import DoubleMLData
 
-from ._utils import _dml_cv_predict, _trimm, _predict_zero_one_propensity, _get_bracket_guess, \
-    _default_kde, _normalize_ipw, _dml_tune, _solve_ipw_score
+from ._utils import (
+    _dml_cv_predict,
+    _trimm,
+    _predict_zero_one_propensity,
+    _get_bracket_guess,
+    _default_kde,
+    _normalize_ipw,
+    _dml_tune,
+    _solve_ipw_score,
+    _cond_targets,
+)
 from ._utils_resampling import DoubleMLResampling
-from ._utils_checks import _check_score, _check_trimming, _check_zero_one_treatment, _check_treatment, \
-    _check_contains_iv, _check_quantile
+from ._utils_checks import (
+    _check_score,
+    _check_trimming,
+    _check_zero_one_treatment,
+    _check_treatment,
+    _check_contains_iv,
+    _check_quantile,
+)
 
 
 class DoubleMLPQ(NonLinearScoreMixin, DoubleML):
@@ -130,21 +144,22 @@ def __init__(self,
             self._kde = _default_kde
         else:
             if not callable(kde):
-                raise TypeError('kde should be either a callable or None. '
-                                '%r was passed.' % kde)
+                raise TypeError("kde should be either a callable or None. "
+                                "%r was passed." % kde)
             self._kde = kde
 
         self._normalize_ipw = normalize_ipw
         self._check_data(self._dml_data)
 
-        valid_score = ['PQ']
+        valid_score = ["PQ"]
         _check_score(self.score, valid_score, allow_callable=False)
         _check_quantile(self.quantile)
         _check_treatment(self.treatment)
 
         if not isinstance(self.normalize_ipw, bool):
-            raise TypeError('Normalization indicator has to be boolean. ' +
-                            f'Object of type {str(type(self.normalize_ipw))} passed.')
+            raise TypeError(
+                "Normalization indicator has to be boolean. " + f"Object of type {str(type(self.normalize_ipw))} passed."
+            )
 
         # initialize starting values and bounds
         self._coef_bounds = (self._dml_data.y.min(), self._dml_data.y.max())
@@ -155,20 +170,23 @@ def __init__(self,
         self._trimming_threshold = trimming_threshold
         _check_trimming(self._trimming_rule, self._trimming_threshold)
 
-        _ = self._check_learner(ml_g, 'ml_g', regressor=False, classifier=True)
-        _ = self._check_learner(ml_m, 'ml_m', regressor=False, classifier=True)
-        self._learner = {'ml_g': ml_g, 'ml_m': ml_m}
-        self._predict_method = {'ml_g': 'predict_proba', 'ml_m': 'predict_proba'}
+        _ = self._check_learner(ml_g, "ml_g", regressor=False, classifier=True)
+        _ = self._check_learner(ml_m, "ml_m", regressor=False, classifier=True)
+        self._learner = {"ml_g": ml_g, "ml_m": ml_m}
+        self._predict_method = {"ml_g": "predict_proba", "ml_m": "predict_proba"}
 
         self._initialize_ml_nuisance_params()
 
         if draw_sample_splitting:
-            obj_dml_resampling = DoubleMLResampling(n_folds=self.n_folds,
-                                                    n_rep=self.n_rep,
-                                                    n_obs=self._dml_data.n_obs,
-                                                    apply_cross_fitting=self.apply_cross_fitting,
-                                                    stratify=self._dml_data.d)
+            obj_dml_resampling = DoubleMLResampling(
+                n_folds=self.n_folds,
+                n_rep=self.n_rep,
+                n_obs=self._dml_data.n_obs,
+                apply_cross_fitting=self.apply_cross_fitting,
+                stratify=self._dml_data.d,
+            )
             self._smpls = obj_dml_resampling.split_samples()
+        self._external_predictions_implemented = True
 
     @property
     def quantile(self):
@@ -214,36 +232,36 @@ def trimming_threshold(self):
 
     @property
     def _score_element_names(self):
-        return ['ind_d', 'g', 'm', 'y']
+        return ["ind_d", "g", "m", "y"]
 
     def _compute_ipw_score(self, theta, d, y, prop):
         score = (d == self.treatment) / prop * (y <= theta) - self.quantile
         return score
 
     def _compute_score(self, psi_elements, coef, inds=None):
-        ind_d = psi_elements['ind_d']
-        g = psi_elements['g']
-        m = psi_elements['m']
-        y = psi_elements['y']
+        ind_d = psi_elements["ind_d"]
+        g = psi_elements["g"]
+        m = psi_elements["m"]
+        y = psi_elements["y"]
 
         if inds is not None:
-            ind_d = psi_elements['ind_d'][inds]
-            g = psi_elements['g'][inds]
-            m = psi_elements['m'][inds]
-            y = psi_elements['y'][inds]
+            ind_d = psi_elements["ind_d"][inds]
+            g = psi_elements["g"][inds]
+            m = psi_elements["m"][inds]
+            y = psi_elements["y"][inds]
 
         score = ind_d * ((y <= coef) - g) / m + g - self.quantile
         return score
 
     def _compute_score_deriv(self, psi_elements, coef, inds=None):
-        ind_d = psi_elements['ind_d']
-        m = psi_elements['m']
-        y = psi_elements['y']
+        ind_d = psi_elements["ind_d"]
+        m = psi_elements["m"]
+        y = psi_elements["y"]
 
         if inds is not None:
-            ind_d = psi_elements['ind_d'][inds]
-            m = psi_elements['m'][inds]
-            y = psi_elements['y'][inds]
+            ind_d = psi_elements["ind_d"][inds]
+            m = psi_elements["m"][inds]
+            y = psi_elements["y"][inds]
 
         score_weights = ind_d / m
 
@@ -253,105 +271,133 @@ def _compute_score_deriv(self, psi_elements, coef, inds=None):
         return deriv
 
     def _initialize_ml_nuisance_params(self):
-        self._params = {learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols}
-                        for learner in ['ml_g', 'ml_m']}
+        self._params = {learner: {key: [None] * self.n_rep for key in self._dml_data.d_cols} for learner in ["ml_g", "ml_m"]}
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
-        x, y = check_X_y(self._dml_data.x, self._dml_data.y,
-                         force_all_finite=False)
-        x, d = check_X_y(x, self._dml_data.d,
-                         force_all_finite=False)
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
+        x, y = check_X_y(self._dml_data.x, self._dml_data.y, force_all_finite=False)
+        x, d = check_X_y(x, self._dml_data.d, force_all_finite=False)
+
+        g_external = external_predictions["ml_g"] is not None
+        m_external = external_predictions["ml_m"] is not None
 
         # initialize nuisance predictions, targets and models
-        g_hat = {'models': None,
-                 'targets': np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
-                 'preds': np.full(shape=self._dml_data.n_obs, fill_value=np.nan)
-                 }
-        m_hat = copy.deepcopy(g_hat)
 
+        if not g_external:
+            g_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
+        if not m_external:
+            m_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+            }
         ipw_vec = np.full(shape=self.n_folds, fill_value=np.nan)
         # initialize models
         fitted_models = {}
         for learner in self.params_names:
             # set nuisance model parameters
-            est_params = self._get_params(learner)
-            if est_params is not None:
-                fitted_models[learner] = [clone(self._learner[learner]).set_params(**est_params[i_fold])
-                                          for i_fold in range(self.n_folds)]
-            else:
-                fitted_models[learner] = [clone(self._learner[learner]) for i_fold in range(self.n_folds)]
+            if (learner == "ml_g" and not g_external) or (learner == "ml_m" and not m_external):
+                est_params = self._get_params(learner)
+                if est_params is not None:
+                    fitted_models[learner] = [
+                        clone(self._learner[learner]).set_params(**est_params[i_fold]) for i_fold in range(self.n_folds)
+                    ]
+                else:
+                    fitted_models[learner] = [clone(self._learner[learner]) for i_fold in range(self.n_folds)]
+        if g_external:
+            g_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_g"],
+            }
+        if m_external:
+            m_hat = {
+                "models": None,
+                "targets": np.full(shape=self._dml_data.n_obs, fill_value=np.nan),
+                "preds": external_predictions["ml_m"],
+            }
 
         # caculate nuisance functions over different folds
-        for i_fold in range(self.n_folds):
-            train_inds = smpls[i_fold][0]
-            test_inds = smpls[i_fold][1]
-
-            # start nested crossfitting
-            train_inds_1, train_inds_2 = train_test_split(train_inds, test_size=0.5,
-                                                          random_state=42, stratify=d[train_inds])
-            smpls_prelim = [(train, test) for train, test in
-                            StratifiedKFold(n_splits=self.n_folds).split(X=train_inds_1, y=d[train_inds_1])]
-
-            d_train_1 = d[train_inds_1]
-            y_train_1 = y[train_inds_1]
-            x_train_1 = x[train_inds_1, :]
-
-            # get a copy of ml_m as a preliminary learner
-            ml_m_prelim = clone(fitted_models['ml_m'][i_fold])
-            m_hat_prelim = _dml_cv_predict(ml_m_prelim, x_train_1, d_train_1,
-                                           method='predict_proba', smpls=smpls_prelim)['preds']
-
-            m_hat_prelim = _trimm(m_hat_prelim, self.trimming_rule, self.trimming_threshold)
-
-            if self._normalize_ipw:
-                m_hat_prelim = _normalize_ipw(m_hat_prelim, d_train_1)
-            if self.treatment == 0:
-                m_hat_prelim = 1 - m_hat_prelim
-
-            # preliminary ipw estimate
-            def ipw_score(theta):
-                res = np.mean(self._compute_ipw_score(theta, d_train_1, y_train_1, m_hat_prelim))
-                return res
-
-            _, bracket_guess = _get_bracket_guess(ipw_score, self._coef_start_val, self._coef_bounds)
-            ipw_est = _solve_ipw_score(ipw_score=ipw_score, bracket_guess=bracket_guess)
-            ipw_vec[i_fold] = ipw_est
-
-            # use the preliminary estimates to fit the nuisance parameters on train_2
-            d_train_2 = d[train_inds_2]
-            y_train_2 = y[train_inds_2]
-            x_train_2 = x[train_inds_2, :]
-
-            dx_treat_train_2 = x_train_2[d_train_2 == self.treatment, :]
-            y_treat_train_2 = y_train_2[d_train_2 == self.treatment]
-
-            fitted_models['ml_g'][i_fold].fit(dx_treat_train_2, y_treat_train_2 <= ipw_est)
-
-            # predict nuisance values on the test data and the corresponding targets
-            g_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_g'][i_fold], x[test_inds, :])
-            g_hat['targets'][test_inds] = y[test_inds] <= ipw_est
-
-            # refit the propensity score on the whole training set
-            fitted_models['ml_m'][i_fold].fit(x[train_inds, :], d[train_inds])
-            m_hat['preds'][test_inds] = _predict_zero_one_propensity(fitted_models['ml_m'][i_fold], x[test_inds, :])
+        if not all([g_external, m_external]):
+            for i_fold in range(self.n_folds):
+                train_inds = smpls[i_fold][0]
+                test_inds = smpls[i_fold][1]
+
+                # start nested crossfitting
+                train_inds_1, train_inds_2 = train_test_split(
+                    train_inds, test_size=0.5, random_state=42, stratify=d[train_inds]
+                )
+                smpls_prelim = [
+                    (train, test)
+                    for train, test in StratifiedKFold(n_splits=self.n_folds).split(X=train_inds_1, y=d[train_inds_1])
+                ]
+
+                d_train_1 = d[train_inds_1]
+                y_train_1 = y[train_inds_1]
+                x_train_1 = x[train_inds_1, :]
+
+                if not m_external:
+                    # get a copy of ml_m as a preliminary learner
+                    ml_m_prelim = clone(fitted_models["ml_m"][i_fold])
+                    m_hat_prelim = _dml_cv_predict(
+                        ml_m_prelim, x_train_1, d_train_1, method="predict_proba", smpls=smpls_prelim
+                    )["preds"]
+                else:
+                    m_hat_prelim = m_hat["preds"][np.concatenate([test for _, test in smpls_prelim])]
+                m_hat_prelim = _trimm(m_hat_prelim, self.trimming_rule, self.trimming_threshold)
+                if self._normalize_ipw:
+                    m_hat_prelim = _normalize_ipw(m_hat_prelim, d_train_1)
+                if self.treatment == 0:
+                    m_hat_prelim = 1 - m_hat_prelim
+
+                # preliminary ipw estimate
+                def ipw_score(theta):
+                    res = np.mean(self._compute_ipw_score(theta, d_train_1, y_train_1, m_hat_prelim))
+                    return res
+
+                _, bracket_guess = _get_bracket_guess(ipw_score, self._coef_start_val, self._coef_bounds)
+                ipw_est = _solve_ipw_score(ipw_score=ipw_score, bracket_guess=bracket_guess)
+                ipw_vec[i_fold] = ipw_est
+
+                # use the preliminary estimates to fit the nuisance parameters on train_2
+                d_train_2 = d[train_inds_2]
+                y_train_2 = y[train_inds_2]
+                x_train_2 = x[train_inds_2, :]
+
+                dx_treat_train_2 = x_train_2[d_train_2 == self.treatment, :]
+                y_treat_train_2 = y_train_2[d_train_2 == self.treatment]
+
+                if not g_external:
+                    fitted_models["ml_g"][i_fold].fit(dx_treat_train_2, y_treat_train_2 <= ipw_est)
+
+                    # predict nuisance values on the test data and the corresponding targets
+                    g_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_g"][i_fold], x[test_inds, :])
+                    g_hat["targets"][test_inds] = y[test_inds] <= ipw_est
+                if not m_external:
+                    # refit the propensity score on the whole training set
+                    fitted_models["ml_m"][i_fold].fit(x[train_inds, :], d[train_inds])
+                    m_hat["preds"][test_inds] = _predict_zero_one_propensity(fitted_models["ml_m"][i_fold], x[test_inds, :])
 
         # set target for propensity score
-        m_hat['targets'] = d
+        m_hat["targets"] = d
 
         # set the target for g to be a float and only relevant values
-        g_hat['targets'] = g_hat['targets'].astype(float)
-        g_hat['targets'][d != self.treatment] = np.nan
+        g_hat["targets"] = _cond_targets(g_hat["targets"], cond_sample=(d == self.treatment))
 
         if return_models:
-            g_hat['models'] = fitted_models['ml_g']
-            m_hat['models'] = fitted_models['ml_m']
+            g_hat["models"] = fitted_models["ml_g"]
+            m_hat["models"] = fitted_models["ml_m"]
 
         # clip propensities and normalize ipw weights
         # this is not done in the score to save computation due to multiple score evaluations
         # to be able to evaluate the raw models the m_hat['preds'] are not changed
-        m_hat_adj = _trimm(m_hat['preds'], self.trimming_rule, self.trimming_threshold)
+
+        m_hat_adj = _trimm(m_hat["preds"], self.trimming_rule, self.trimming_threshold)
         if self._normalize_ipw:
-            if self.dml_procedure == 'dml1':
+            if self.dml_procedure == "dml1":
                 for _, test_index in smpls:
                     m_hat_adj[test_index] = _normalize_ipw(m_hat_adj[test_index], d[test_index])
             else:
@@ -359,63 +405,74 @@ def ipw_score(theta):
 
         if self.treatment == 0:
             m_hat_adj = 1 - m_hat_adj
-
         # readjust start value for minimization
-        self._coef_start_val = np.mean(ipw_vec)
-
-        psi_elements = {'ind_d': d == self.treatment, 'g': g_hat['preds'],
-                        'm': m_hat_adj, 'y': y}
-
-        preds = {'predictions': {'ml_g': g_hat['preds'],
-                                 'ml_m': m_hat['preds']},
-                 'targets': {'ml_g': g_hat['targets'],
-                             'ml_m': m_hat['targets']},
-                 'models': {'ml_g': g_hat['models'],
-                            'ml_m': m_hat['models']}
-                 }
+        if not g_external or not m_external:
+            self._coef_start_val = np.mean(ipw_vec)
+
+        psi_elements = {"ind_d": d == self.treatment, "g": g_hat["preds"], "m": m_hat_adj, "y": y}
+
+        preds = {
+            "predictions": {"ml_g": g_hat["preds"], "ml_m": m_hat["preds"]},
+            "targets": {"ml_g": g_hat["targets"], "ml_m": m_hat["targets"]},
+            "models": {"ml_g": g_hat["models"], "ml_m": m_hat["models"]},
+        }
         return psi_elements, preds
 
-    def _nuisance_tuning(self, smpls, param_grids, scoring_methods, n_folds_tune, n_jobs_cv,
-                         search_mode, n_iter_randomized_search):
-        x, y = check_X_y(self._dml_data.x, self._dml_data.y,
-                         force_all_finite=False)
-        x, d = check_X_y(x, self._dml_data.d,
-                         force_all_finite=False)
+    def _nuisance_tuning(
+        self, smpls, param_grids, scoring_methods, n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search
+    ):
+        x, y = check_X_y(self._dml_data.x, self._dml_data.y, force_all_finite=False)
+        x, d = check_X_y(x, self._dml_data.d, force_all_finite=False)
 
         if scoring_methods is None:
-            scoring_methods = {'ml_g': None,
-                               'ml_m': None}
+            scoring_methods = {"ml_g": None, "ml_m": None}
 
         train_inds = [train_index for (train_index, _) in smpls]
         train_inds_treat = [np.intersect1d(np.where(d == self.treatment)[0], train) for train, _ in smpls]
 
         # use self._coef_start_val as a very crude approximation of ipw_est
         approx_goal = y <= np.quantile(y[d == self.treatment], self.quantile)
-        g_tune_res = _dml_tune(approx_goal, x, train_inds_treat,
-                               self._learner['ml_g'], param_grids['ml_g'], scoring_methods['ml_g'],
-                               n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
-
-        m_tune_res = _dml_tune(d, x, train_inds,
-                               self._learner['ml_m'], param_grids['ml_m'], scoring_methods['ml_m'],
-                               n_folds_tune, n_jobs_cv, search_mode, n_iter_randomized_search)
+        g_tune_res = _dml_tune(
+            approx_goal,
+            x,
+            train_inds_treat,
+            self._learner["ml_g"],
+            param_grids["ml_g"],
+            scoring_methods["ml_g"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
+
+        m_tune_res = _dml_tune(
+            d,
+            x,
+            train_inds,
+            self._learner["ml_m"],
+            param_grids["ml_m"],
+            scoring_methods["ml_m"],
+            n_folds_tune,
+            n_jobs_cv,
+            search_mode,
+            n_iter_randomized_search,
+        )
 
         g_best_params = [xx.best_params_ for xx in g_tune_res]
         m_best_params = [xx.best_params_ for xx in m_tune_res]
 
-        params = {'ml_g': g_best_params,
-                  'ml_m': m_best_params}
-        tune_res = {'g_tune': g_tune_res,
-                    'm_tune': m_tune_res}
+        params = {"ml_g": g_best_params, "ml_m": m_best_params}
+        tune_res = {"g_tune": g_tune_res, "m_tune": m_tune_res}
 
-        res = {'params': params,
-               'tune_res': tune_res}
+        res = {"params": params, "tune_res": tune_res}
 
         return res
 
     def _check_data(self, obj_dml_data):
         if not isinstance(obj_dml_data, DoubleMLData):
-            raise TypeError('The data must be of DoubleMLData type. '
-                            f'{str(obj_dml_data)} of type {str(type(obj_dml_data))} was passed.')
+            raise TypeError(
+                "The data must be of DoubleMLData type. " f"{str(obj_dml_data)} of type {str(type(obj_dml_data))} was passed."
+            )
         _check_contains_iv(obj_dml_data)
         _check_zero_one_treatment(self)
         return
diff --git a/doubleml/double_ml_qte.py b/doubleml/double_ml_qte.py
index c640abf7..9633434f 100644
--- a/doubleml/double_ml_qte.py
+++ b/doubleml/double_ml_qte.py
@@ -161,7 +161,7 @@ def __init__(self,
         self._modellist_0, self._modellist_1 = self._initialize_models()
 
         # initialize arrays according to obj_dml_data and the resampling settings
-        self._psi0, self._psi1, self._psi0_deriv, self._psi1_deriv,\
+        self._psi0, self._psi1, self._psi0_deriv, self._psi1_deriv, \
             self._coef, self._se, self._all_coef, self._all_se, self._all_dml1_coef = self._initialize_arrays()
 
         # also initialize bootstrap arrays with the default number of bootstrap replications
@@ -386,7 +386,7 @@ def __psi1_deriv(self):
     def __all_se(self):
         return self._all_se[self._i_quant, self._i_rep]
 
-    def fit(self, n_jobs_models=None, n_jobs_cv=None, store_predictions=True, store_models=False):
+    def fit(self, n_jobs_models=None, n_jobs_cv=None, store_predictions=True, store_models=False, external_predictions=None):
         """
         Estimate DoubleMLQTE models.
 
@@ -415,6 +415,9 @@ def fit(self, n_jobs_models=None, n_jobs_cv=None, store_predictions=True, store_
         self : object
         """
 
+        if external_predictions is not None:
+            raise NotImplementedError(f"External predictions not implemented for {self.__class__.__name__}.")
+
         # parallel estimation of the quantiles
         parallel = Parallel(n_jobs=n_jobs_models, verbose=0, pre_dispatch='2*n_jobs')
         fitted_models = parallel(delayed(self._fit_quantile)(i_quant, n_jobs_cv, store_predictions, store_models)
diff --git a/doubleml/utils/__init__.py b/doubleml/utils/__init__.py
new file mode 100644
index 00000000..b3fbb9f0
--- /dev/null
+++ b/doubleml/utils/__init__.py
@@ -0,0 +1,7 @@
+from .dummy_learners import dummy_classifier
+from .dummy_learners import dummy_regressor
+
+__all__ = [
+    "dummy_classifier",
+    "dummy_regressor",
+]
diff --git a/doubleml/utils/dummy_learners.py b/doubleml/utils/dummy_learners.py
new file mode 100644
index 00000000..4d771b20
--- /dev/null
+++ b/doubleml/utils/dummy_learners.py
@@ -0,0 +1,66 @@
+from sklearn.base import BaseEstimator
+
+
+class dummy_regressor(BaseEstimator):
+    """
+    A dummy regressor that raises an AttributeError when attempting to access
+    its fit, predict, or set_params methods.
+    Attributes
+    ----------
+    _estimator_type : str
+        Type of the estimator, set to "regressor".
+    Methods
+    -------
+    fit(*args)
+        Raises AttributeError: "Accessed fit method of DummyRegressor!"
+    predict(*args)
+        Raises AttributeError: "Accessed predict method of DummyRegressor!"
+    set_params(*args)
+        Raises AttributeError: "Accessed set_params method of DummyRegressor!"
+    """
+
+    _estimator_type = "regressor"
+
+    def fit(*args):
+        raise AttributeError("Accessed fit method of dummy_regressor!")
+
+    def predict(*args):
+        raise AttributeError("Accessed predict method of dummy_regressor!")
+
+    def set_params(*args):
+        raise AttributeError("Accessed set_params method of dummy_regressor!")
+
+
+class dummy_classifier(BaseEstimator):
+    """
+    A dummy classifier that raises an AttributeError when attempting to access
+    its fit, predict, set_params, or predict_proba methods.
+    Attributes
+    ----------
+    _estimator_type : str
+        Type of the estimator, set to "classifier".
+    Methods
+    -------
+    fit(*args)
+        Raises AttributeError: "Accessed fit method of DummyClassifier!"
+    predict(*args)
+        Raises AttributeError: "Accessed predict method of DummyClassifier!"
+    set_params(*args)
+        Raises AttributeError: "Accessed set_params method of DummyClassifier!"
+    predict_proba(*args, **kwargs)
+        Raises AttributeError: "Accessed predict_proba method of DummyClassifier!"
+    """
+
+    _estimator_type = "classifier"
+
+    def fit(*args):
+        raise AttributeError("Accessed fit method of dummy_classifier!")
+
+    def predict(*args):
+        raise AttributeError("Accessed predict method of dummy_classifier!")
+
+    def set_params(*args):
+        raise AttributeError("Accessed set_params method of dummy_classifier!")
+
+    def predict_proba(*args, **kwargs):
+        raise AttributeError("Accessed predict_proba method of dummy_classifier!")


Test patch:
diff --git a/doubleml/tests/_utils_irm_manual.py b/doubleml/tests/_utils_irm_manual.py
index 5f162bee..c46fe605 100644
--- a/doubleml/tests/_utils_irm_manual.py
+++ b/doubleml/tests/_utils_irm_manual.py
@@ -130,20 +130,23 @@ def irm_dml1(y, x, d, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls, s
     u_hat0, u_hat1, g_hat0, g_hat1, m_hat, p_hat = compute_iivm_residuals(
         y, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls)
 
+    m_hat_adj = np.full_like(m_hat, np.nan, dtype='float64')
     if normalize_ipw:
         for _, test_index in smpls:
-            m_hat[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
+            m_hat_adj[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
+    else:
+        m_hat_adj = m_hat
 
     for idx, (_, test_index) in enumerate(smpls):
         thetas[idx] = irm_orth(g_hat0[test_index], g_hat1[test_index],
-                               m_hat[test_index], p_hat[test_index],
+                               m_hat_adj[test_index], p_hat[test_index],
                                u_hat0[test_index], u_hat1[test_index],
                                d[test_index], score)
     theta_hat = np.mean(thetas)
 
     if len(smpls) > 1:
         se = np.sqrt(var_irm(theta_hat, g_hat0, g_hat1,
-                             m_hat, p_hat,
+                             m_hat_adj, p_hat,
                              u_hat0, u_hat1,
                              d, score, n_obs))
     else:
@@ -151,7 +154,7 @@ def irm_dml1(y, x, d, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls, s
         test_index = smpls[0][1]
         n_obs = len(test_index)
         se = np.sqrt(var_irm(theta_hat, g_hat0[test_index], g_hat1[test_index],
-                             m_hat[test_index], p_hat[test_index],
+                             m_hat_adj[test_index], p_hat[test_index],
                              u_hat0[test_index], u_hat1[test_index],
                              d[test_index], score, n_obs))
 
@@ -164,12 +167,14 @@ def irm_dml2(y, x, d, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls, s
         y, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls)
 
     if normalize_ipw:
-        m_hat = _normalize_ipw(m_hat, d)
+        m_hat_adj = _normalize_ipw(m_hat, d)
+    else:
+        m_hat_adj = m_hat
 
-    theta_hat = irm_orth(g_hat0, g_hat1, m_hat, p_hat,
+    theta_hat = irm_orth(g_hat0, g_hat1, m_hat_adj, p_hat,
                          u_hat0, u_hat1, d, score)
     se = np.sqrt(var_irm(theta_hat, g_hat0, g_hat1,
-                         m_hat, p_hat,
+                         m_hat_adj, p_hat,
                          u_hat0, u_hat1,
                          d, score, n_obs))
 
@@ -240,12 +245,15 @@ def boot_irm_single_split(theta, y, d, g_hat0_list, g_hat1_list, m_hat_list, p_h
     u_hat0, u_hat1, g_hat0, g_hat1, m_hat, p_hat = compute_iivm_residuals(
         y, g_hat0_list, g_hat1_list, m_hat_list, p_hat_list, smpls)
 
+    m_hat_adj = np.full_like(m_hat, np.nan, dtype='float64')
     if normalize_ipw:
         if dml_procedure == 'dml1':
             for _, test_index in smpls:
-                m_hat[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
+                m_hat_adj[test_index] = _normalize_ipw(m_hat[test_index], d[test_index])
         else:
-            m_hat = _normalize_ipw(m_hat, d)
+            m_hat_adj = _normalize_ipw(m_hat, d)
+    else:
+        m_hat_adj = m_hat
 
     if apply_cross_fitting:
         if score == 'ATE':
@@ -263,13 +271,13 @@ def boot_irm_single_split(theta, y, d, g_hat0_list, g_hat1_list, m_hat_list, p_h
 
     if score == 'ATE':
         psi = g_hat1 - g_hat0 \
-                + np.divide(np.multiply(d, u_hat1), m_hat) \
-                - np.divide(np.multiply(1.-d, u_hat0), 1.-m_hat) - theta
+                + np.divide(np.multiply(d, u_hat1), m_hat_adj) \
+                - np.divide(np.multiply(1.-d, u_hat0), 1.-m_hat_adj) - theta
     else:
         assert score == 'ATTE'
         psi = np.divide(np.multiply(d, u_hat0), p_hat) \
-            - np.divide(np.multiply(m_hat, np.multiply(1.-d, u_hat0)),
-                        np.multiply(p_hat, (1.-m_hat))) \
+            - np.divide(np.multiply(m_hat_adj, np.multiply(1.-d, u_hat0)),
+                        np.multiply(p_hat, (1.-m_hat_adj))) \
             - theta * np.divide(d, p_hat)
 
     boot_theta, boot_t_stat = boot_manual(psi, J, smpls, se, weights, n_rep_boot, apply_cross_fitting)
diff --git a/doubleml/tests/test_did_external_predictions.py b/doubleml/tests/test_did_external_predictions.py
new file mode 100644
index 00000000..43138cf8
--- /dev/null
+++ b/doubleml/tests/test_did_external_predictions.py
@@ -0,0 +1,59 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression, LogisticRegression
+from doubleml import DoubleMLDID
+from doubleml.datasets import make_did_SZ2020
+from doubleml.utils import dummy_regressor, dummy_classifier
+from ._utils import draw_smpls
+
+
+@pytest.fixture(scope="module", params=["observational", "experimental"])
+def did_score(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_did_fixture(did_score, dml_procedure, n_rep):
+    ext_predictions = {"d": {}}
+    dml_data = make_did_SZ2020(n_obs=500, return_type="DoubleMLData")
+    all_smpls = draw_smpls(len(dml_data.y), 5, n_rep=n_rep, groups=dml_data.d)
+    kwargs = {
+        "obj_dml_data": dml_data,
+        "score": did_score,
+        "n_rep": n_rep,
+        "dml_procedure": dml_procedure,
+        "draw_sample_splitting": False
+    }
+    dml_did = DoubleMLDID(ml_g=LinearRegression(), ml_m=LogisticRegression(), **kwargs)
+    dml_did.set_sample_splitting(all_smpls)
+    np.random.seed(3141)
+    dml_did.fit(store_predictions=True)
+
+    ext_predictions["d"]["ml_g0"] = dml_did.predictions["ml_g0"][:, :, 0]
+    ext_predictions["d"]["ml_g1"] = dml_did.predictions["ml_g1"][:, :, 0]
+    ext_predictions["d"]["ml_m"] = dml_did.predictions["ml_m"][:, :, 0]
+
+    dml_did_ext = DoubleMLDID(ml_g=dummy_regressor(), ml_m=dummy_classifier(), **kwargs)
+    dml_did_ext.set_sample_splitting(all_smpls)
+    np.random.seed(3141)
+    dml_did_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_did.coef[0], "coef_ext": dml_did_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_did_coef(doubleml_did_fixture):
+    assert math.isclose(doubleml_did_fixture["coef_normal"], doubleml_did_fixture["coef_ext"], rel_tol=1e-9, abs_tol=1e-3)
diff --git a/doubleml/tests/test_didcs_external_predictions.py b/doubleml/tests/test_didcs_external_predictions.py
new file mode 100644
index 00000000..a1ffda5e
--- /dev/null
+++ b/doubleml/tests/test_didcs_external_predictions.py
@@ -0,0 +1,62 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression, LogisticRegression
+from doubleml import DoubleMLDIDCS
+from doubleml.datasets import make_did_SZ2020
+from doubleml.utils import dummy_regressor, dummy_classifier
+from ._utils import draw_smpls
+
+
+@pytest.fixture(scope="module", params=["observational", "experimental"])
+def did_score(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_didcs_fixture(did_score, dml_procedure, n_rep):
+    ext_predictions = {"d": {}}
+    dml_data = make_did_SZ2020(n_obs=500, cross_sectional_data=True, return_type="DoubleMLData")
+    all_smpls = draw_smpls(len(dml_data.y), 5, n_rep=n_rep, groups=dml_data.d)
+    kwargs = {
+        "obj_dml_data": dml_data,
+        "score": did_score,
+        "n_rep": n_rep,
+        "n_folds": 5,
+        "dml_procedure": dml_procedure,
+        "draw_sample_splitting": False
+    }
+    dml_did_cs = DoubleMLDIDCS(ml_g=LinearRegression(), ml_m=LogisticRegression(), **kwargs)
+    dml_did_cs.set_sample_splitting(all_smpls)
+    np.random.seed(3141)
+    dml_did_cs.fit(store_predictions=True)
+
+    ext_predictions["d"]["ml_g_d0_t0"] = dml_did_cs.predictions["ml_g_d0_t0"][:, :, 0]
+    ext_predictions["d"]["ml_g_d0_t1"] = dml_did_cs.predictions["ml_g_d0_t1"][:, :, 0]
+    ext_predictions["d"]["ml_g_d1_t0"] = dml_did_cs.predictions["ml_g_d1_t0"][:, :, 0]
+    ext_predictions["d"]["ml_g_d1_t1"] = dml_did_cs.predictions["ml_g_d1_t1"][:, :, 0]
+    ext_predictions["d"]["ml_m"] = dml_did_cs.predictions["ml_m"][:, :, 0]
+
+    dml_did_cs_ext = DoubleMLDIDCS(ml_g=dummy_regressor(), ml_m=dummy_classifier(), **kwargs)
+    dml_did_cs_ext.set_sample_splitting(all_smpls)
+    np.random.seed(3141)
+    dml_did_cs_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_did_cs.coef[0], "coef_ext": dml_did_cs_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_didcs_coef(doubleml_didcs_fixture):
+    assert math.isclose(doubleml_didcs_fixture["coef_normal"], doubleml_didcs_fixture["coef_ext"], rel_tol=1e-9, abs_tol=1e-3)
diff --git a/doubleml/tests/test_doubleml_exceptions.py b/doubleml/tests/test_doubleml_exceptions.py
index bb5a75bd..b840cc9b 100644
--- a/doubleml/tests/test_doubleml_exceptions.py
+++ b/doubleml/tests/test_doubleml_exceptions.py
@@ -1402,3 +1402,84 @@ def test_doubleml_exception_policytree():
     msg = 'Only implemented for one repetition. Number of repetitions is 2.'
     with pytest.raises(NotImplementedError, match=msg):
         dml_irm_obj.policy_tree(features=2, depth=1)
+
+
+@pytest.mark.ci
+def test_double_ml_external_predictions():
+    dml_irm_obj = DoubleMLIRM(dml_data_irm,
+                              ml_g=Lasso(),
+                              ml_m=LogisticRegression(),
+                              trimming_threshold=0.05,
+                              n_folds=5,
+                              score='ATE',
+                              n_rep=2)
+
+    msg = "external_predictions must be a dictionary. ml_m of type <class 'str'> was passed."
+    with pytest.raises(TypeError, match=msg):
+        dml_irm_obj.fit(external_predictions="ml_m")
+
+    dml_irm_obj = DoubleMLIRM(dml_data_irm,
+                              ml_g=Lasso(),
+                              ml_m=LogisticRegression(),
+                              trimming_threshold=0.05,
+                              n_folds=5,
+                              score='ATE',
+                              n_rep=1)
+
+    predictions = {'d': 'test', 'd_f': 'test'}
+    msg = (r"Invalid external_predictions. Invalid treatment variable in \['d', 'd_f'\]. "
+           "Valid treatment variables d.")
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': 'test'}
+    msg = ("external_predictions must be a nested dictionary. "
+           "For treatment d a value of type <class 'str'> was passed.")
+    with pytest.raises(TypeError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_f': 'test'}}
+    msg = ("Invalid external_predictions. "
+           r"Invalid nuisance learner for treatment d in \['ml_f'\]. "
+           "Valid nuisance learners ml_g0 or ml_g1 or ml_m.")
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_m': 'test', 'ml_f': 'test'}}
+    msg = ("Invalid external_predictions. "
+           r"Invalid nuisance learner for treatment d in \['ml_m', 'ml_f'\]. "
+           "Valid nuisance learners ml_g0 or ml_g1 or ml_m.")
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_m': 'test'}}
+    msg = ("Invalid external_predictions. "
+           "The values of the nested list must be a numpy array. "
+           "Invalid predictions for treatment d and learner ml_m. "
+           "Object of type <class 'str'> was passed.")
+    with pytest.raises(TypeError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_m': np.array([0])}}
+    msg = ('Invalid external_predictions. '
+           r'The supplied predictions have to be of shape \(100, 1\). '
+           'Invalid predictions for treatment d and learner ml_m. '
+           r'Predictions of shape \(1,\) passed.')
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_m': np.zeros(100)}}
+    msg = ('Invalid external_predictions. '
+           r'The supplied predictions have to be of shape \(100, 1\). '
+           'Invalid predictions for treatment d and learner ml_m. '
+           r'Predictions of shape \(100,\) passed.')
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
+
+    predictions = {'d': {'ml_m': np.ones(shape=(5, 3))}}
+    msg = ('Invalid external_predictions. '
+           r'The supplied predictions have to be of shape \(100, 1\). '
+           'Invalid predictions for treatment d and learner ml_m. '
+           r'Predictions of shape \(5, 3\) passed.')
+    with pytest.raises(ValueError, match=msg):
+        dml_irm_obj.fit(external_predictions=predictions)
diff --git a/doubleml/tests/test_doubleml_exceptions_ext_preds.py b/doubleml/tests/test_doubleml_exceptions_ext_preds.py
new file mode 100644
index 00000000..4be4430f
--- /dev/null
+++ b/doubleml/tests/test_doubleml_exceptions_ext_preds.py
@@ -0,0 +1,23 @@
+import pytest
+from doubleml import DoubleMLCVAR, DoubleMLQTE, DoubleMLData
+from doubleml.datasets import make_irm_data
+from doubleml.utils import dummy_regressor, dummy_classifier
+
+df_irm = make_irm_data(n_obs=10, dim_x=2, theta=0.5, return_type="DataFrame")
+ext_predictions = {"d": {}}
+
+
+@pytest.mark.ci
+def test_cvar_external_prediction_exception():
+    msg = "External predictions not implemented for DoubleMLCVAR."
+    with pytest.raises(NotImplementedError, match=msg):
+        cvar = DoubleMLCVAR(DoubleMLData(df_irm, "y", "d"), dummy_regressor(), dummy_classifier(), treatment=1)
+        cvar.fit(external_predictions=ext_predictions)
+
+
+@pytest.mark.ci
+def test_qte_external_prediction_exception():
+    msg = "External predictions not implemented for DoubleMLQTE."
+    with pytest.raises(NotImplementedError, match=msg):
+        qte = DoubleMLQTE(DoubleMLData(df_irm, "y", "d"), dummy_classifier(), dummy_classifier())
+        qte.fit(external_predictions=ext_predictions)
diff --git a/doubleml/tests/test_dummy_learners.py b/doubleml/tests/test_dummy_learners.py
new file mode 100644
index 00000000..a357345c
--- /dev/null
+++ b/doubleml/tests/test_dummy_learners.py
@@ -0,0 +1,46 @@
+import pytest
+import numpy as np
+from doubleml.utils import dummy_regressor, dummy_classifier
+from sklearn.base import clone
+
+
+@pytest.fixture(scope="module")
+def dl_fixture():
+    fixture = {
+        "dummy_regressor": dummy_regressor(),
+        "dummy_classifier": dummy_classifier(),
+        "X": np.random.normal(0, 1, size=(100, 10)),
+        "y_con": np.random.normal(0, 1, size=(100, 1)),
+        "y_cat": np.random.binomial(1, 0.5, size=(100, 1)),
+    }
+
+    return fixture
+
+
+@pytest.mark.ci
+def test_fit(dl_fixture):
+    msg = "Accessed fit method of dummy_regressor!"
+    with pytest.raises(AttributeError, match=msg):
+        dl_fixture["dummy_regressor"].fit(dl_fixture["X"], dl_fixture["y_con"])
+    msg = "Accessed fit method of dummy_classifier!"
+    with pytest.raises(AttributeError, match=msg):
+        dl_fixture["dummy_classifier"].fit(dl_fixture["X"], dl_fixture["y_cat"])
+
+
+@pytest.mark.ci
+def test_predict(dl_fixture):
+    msg = "Accessed predict method of dummy_regressor!"
+    with pytest.raises(AttributeError, match=msg):
+        dl_fixture["dummy_regressor"].predict(dl_fixture["X"])
+    msg = "Accessed predict method of dummy_classifier!"
+    with pytest.raises(AttributeError, match=msg):
+        dl_fixture["dummy_classifier"].predict(dl_fixture["X"])
+
+
+@pytest.mark.ci
+def test_clone(dl_fixture):
+    try:
+        _ = clone(dl_fixture["dummy_regressor"])
+        _ = clone(dl_fixture["dummy_classifier"])
+    except Exception as e:
+        pytest.fail(f"clone() raised an exception:\n{str(e)}\n")
diff --git a/doubleml/tests/test_iivm_external_predictions.py b/doubleml/tests/test_iivm_external_predictions.py
new file mode 100644
index 00000000..548f2297
--- /dev/null
+++ b/doubleml/tests/test_iivm_external_predictions.py
@@ -0,0 +1,74 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression, LogisticRegression
+from doubleml import DoubleMLIIVM, DoubleMLData
+from doubleml.datasets import make_iivm_data
+from doubleml.utils import dummy_regressor, dummy_classifier
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def adapted_doubleml_fixture(dml_procedure, n_rep):
+    ext_predictions = {"d": {}}
+
+    data = make_iivm_data(
+        n_obs=500, dim_x=20, theta=0.5, alpha_x=1.0, return_type="DataFrame"
+    )
+
+    np.random.seed(3141)
+
+    dml_data = DoubleMLData(data, "y", "d", z_cols="z")
+
+    kwargs = {
+        "obj_dml_data": dml_data,
+        "score": "LATE",
+        "n_rep": n_rep,
+        "dml_procedure": dml_procedure,
+    }
+
+    dml_iivm = DoubleMLIIVM(
+        ml_g=LinearRegression(),
+        ml_m=LogisticRegression(),
+        ml_r=LogisticRegression(),
+        **kwargs,
+    )
+    np.random.seed(3141)
+
+    dml_iivm.fit(store_predictions=True)
+
+    ext_predictions["d"]["ml_g0"] = dml_iivm.predictions["ml_g0"][:, :, 0]
+    ext_predictions["d"]["ml_g1"] = dml_iivm.predictions["ml_g1"][:, :, 0]
+    ext_predictions["d"]["ml_m"] = dml_iivm.predictions["ml_m"][:, :, 0]
+    ext_predictions["d"]["ml_r0"] = dml_iivm.predictions["ml_r0"][:, :, 0]
+    ext_predictions["d"]["ml_r1"] = dml_iivm.predictions["ml_r1"][:, :, 0]
+
+    dml_iivm_ext = DoubleMLIIVM(
+        ml_g=dummy_regressor(), ml_m=dummy_classifier(), ml_r=dummy_classifier(), **kwargs
+    )
+
+    np.random.seed(3141)
+    dml_iivm_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_iivm.coef[0], "coef_ext": dml_iivm_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_adapted_doubleml_coef(adapted_doubleml_fixture):
+    assert math.isclose(
+        adapted_doubleml_fixture["coef_normal"],
+        adapted_doubleml_fixture["coef_ext"],
+        rel_tol=1e-9,
+        abs_tol=1e-4,
+    )
diff --git a/doubleml/tests/test_irm.py b/doubleml/tests/test_irm.py
index e4842ffa..73dd5cca 100644
--- a/doubleml/tests/test_irm.py
+++ b/doubleml/tests/test_irm.py
@@ -87,10 +87,31 @@ def dml_irm_fixture(generate_data_irm, learner, score, dml_procedure, normalize_
                          normalize_ipw=normalize_ipw,
                          trimming_threshold=trimming_threshold)
 
+    np.random.seed(3141)
+    # test with external nuisance predictions
+    dml_irm_obj_ext = dml.DoubleMLIRM(obj_dml_data,
+                                      ml_g, ml_m,
+                                      n_folds,
+                                      score=score,
+                                      dml_procedure=dml_procedure,
+                                      normalize_ipw=normalize_ipw,
+                                      draw_sample_splitting=False,
+                                      trimming_threshold=trimming_threshold)
+
+    # synchronize the sample splitting
+    dml_irm_obj_ext.set_sample_splitting(all_smpls=all_smpls)
+
+    prediction_dict = {'d': {'ml_g0': dml_irm_obj.predictions['ml_g0'].reshape(-1, 1),
+                             'ml_g1': dml_irm_obj.predictions['ml_g1'].reshape(-1, 1),
+                             'ml_m': dml_irm_obj.predictions['ml_m'].reshape(-1, 1)}}
+    dml_irm_obj_ext.fit(external_predictions=prediction_dict)
+
     res_dict = {'coef': dml_irm_obj.coef,
                 'coef_manual': res_manual['theta'],
+                'coef_ext': dml_irm_obj_ext.coef,
                 'se': dml_irm_obj.se,
                 'se_manual': res_manual['se'],
+                'se_ext': dml_irm_obj_ext.se,
                 'boot_methods': boot_methods}
 
     for bootstrap in boot_methods:
@@ -104,10 +125,14 @@ def dml_irm_fixture(generate_data_irm, learner, score, dml_procedure, normalize_
 
         np.random.seed(3141)
         dml_irm_obj.bootstrap(method=bootstrap, n_rep_boot=n_rep_boot)
+        np.random.seed(3141)
+        dml_irm_obj_ext.bootstrap(method=bootstrap, n_rep_boot=n_rep_boot)
         res_dict['boot_coef' + bootstrap] = dml_irm_obj.boot_coef
         res_dict['boot_t_stat' + bootstrap] = dml_irm_obj.boot_t_stat
         res_dict['boot_coef' + bootstrap + '_manual'] = boot_theta
         res_dict['boot_t_stat' + bootstrap + '_manual'] = boot_t_stat
+        res_dict['boot_coef' + bootstrap + '_ext'] = dml_irm_obj_ext.boot_coef
+        res_dict['boot_t_stat' + bootstrap + '_ext'] = dml_irm_obj_ext.boot_t_stat
 
     # sensitivity tests
     res_dict['sensitivity_elements'] = dml_irm_obj.sensitivity_elements
@@ -125,16 +150,22 @@ def dml_irm_fixture(generate_data_irm, learner, score, dml_procedure, normalize_
 
 @pytest.mark.ci
 def test_dml_irm_coef(dml_irm_fixture):
-    assert math.isclose(dml_irm_fixture['coef'],
+    assert math.isclose(dml_irm_fixture['coef'][0],
                         dml_irm_fixture['coef_manual'],
                         rel_tol=1e-9, abs_tol=1e-4)
+    assert math.isclose(dml_irm_fixture['coef'][0],
+                        dml_irm_fixture['coef_ext'][0],
+                        rel_tol=1e-9, abs_tol=1e-4)
 
 
 @pytest.mark.ci
 def test_dml_irm_se(dml_irm_fixture):
-    assert math.isclose(dml_irm_fixture['se'],
+    assert math.isclose(dml_irm_fixture['se'][0],
                         dml_irm_fixture['se_manual'],
                         rel_tol=1e-9, abs_tol=1e-4)
+    assert math.isclose(dml_irm_fixture['se'][0],
+                        dml_irm_fixture['se_ext'][0],
+                        rel_tol=1e-9, abs_tol=1e-4)
 
 
 @pytest.mark.ci
@@ -143,9 +174,15 @@ def test_dml_irm_boot(dml_irm_fixture):
         assert np.allclose(dml_irm_fixture['boot_coef' + bootstrap],
                            dml_irm_fixture['boot_coef' + bootstrap + '_manual'],
                            rtol=1e-9, atol=1e-4)
+        assert np.allclose(dml_irm_fixture['boot_coef' + bootstrap],
+                           dml_irm_fixture['boot_coef' + bootstrap + '_ext'],
+                           rtol=1e-9, atol=1e-4)
         assert np.allclose(dml_irm_fixture['boot_t_stat' + bootstrap],
                            dml_irm_fixture['boot_t_stat' + bootstrap + '_manual'],
                            rtol=1e-9, atol=1e-4)
+        assert np.allclose(dml_irm_fixture['boot_t_stat' + bootstrap],
+                           dml_irm_fixture['boot_t_stat' + bootstrap + '_ext'],
+                           rtol=1e-9, atol=1e-4)
 
 
 @pytest.mark.ci
diff --git a/doubleml/tests/test_irm_external_predictions.py b/doubleml/tests/test_irm_external_predictions.py
new file mode 100644
index 00000000..533e6250
--- /dev/null
+++ b/doubleml/tests/test_irm_external_predictions.py
@@ -0,0 +1,77 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression, LogisticRegression
+from doubleml import DoubleMLIRM, DoubleMLData
+from doubleml.datasets import make_irm_data
+from doubleml.utils import dummy_regressor, dummy_classifier
+
+
+@pytest.fixture(scope="module", params=["ATE", "ATTE"])
+def irm_score(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_m_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_g_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_irm_fixture(irm_score, dml_procedure, n_rep, set_ml_m_ext, set_ml_g_ext):
+    ext_predictions = {"d": {}}
+
+    x, y, d = make_irm_data(n_obs=500, dim_x=20, theta=0.5, return_type="np.array")
+
+    np.random.seed(3141)
+
+    dml_data = DoubleMLData.from_arrays(x=x, y=y, d=d)
+
+    kwargs = {"obj_dml_data": dml_data, "score": irm_score, "n_rep": n_rep, "dml_procedure": dml_procedure}
+
+    dml_irm = DoubleMLIRM(ml_g=LinearRegression(), ml_m=LogisticRegression(), **kwargs)
+    np.random.seed(3141)
+
+    dml_irm.fit(store_predictions=True)
+
+    if set_ml_m_ext:
+        ext_predictions["d"]["ml_m"] = dml_irm.predictions["ml_m"][:, :, 0]
+        ml_m = dummy_classifier()
+    else:
+        ml_m = LogisticRegression(random_state=42)
+
+    if set_ml_g_ext:
+        ext_predictions["d"]["ml_g0"] = dml_irm.predictions["ml_g0"][:, :, 0]
+        ext_predictions["d"]["ml_g1"] = dml_irm.predictions["ml_g1"][:, :, 0]
+        ml_g = dummy_regressor()
+    else:
+        ml_g = LinearRegression()
+
+    dml_irm_ext = DoubleMLIRM(ml_g=ml_g, ml_m=ml_m, **kwargs)
+
+    np.random.seed(3141)
+    dml_irm_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_irm.coef[0], "coef_ext": dml_irm_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_irm_coef(doubleml_irm_fixture):
+    assert math.isclose(doubleml_irm_fixture["coef_normal"], doubleml_irm_fixture["coef_ext"], rel_tol=1e-9, abs_tol=1e-4)
diff --git a/doubleml/tests/test_lpq.py b/doubleml/tests/test_lpq.py
index e7550e06..beb8b6a0 100644
--- a/doubleml/tests/test_lpq.py
+++ b/doubleml/tests/test_lpq.py
@@ -119,9 +119,9 @@ def dml_lpq_fixture(generate_data_local_quantiles, treatment, quantile, learner,
                              normalize_ipw=normalize_ipw, kde=kde,
                              n_rep=1, trimming_threshold=trimming_threshold)
 
-    res_dict = {'coef': dml_lpq_obj.coef,
+    res_dict = {'coef': dml_lpq_obj.coef[0],
                 'coef_manual': res_manual['lpq'],
-                'se': dml_lpq_obj.se,
+                'se': dml_lpq_obj.se[0],
                 'se_manual': res_manual['se']}
 
     return res_dict
diff --git a/doubleml/tests/test_lpq_external_predictions.py b/doubleml/tests/test_lpq_external_predictions.py
new file mode 100644
index 00000000..af46191f
--- /dev/null
+++ b/doubleml/tests/test_lpq_external_predictions.py
@@ -0,0 +1,72 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LogisticRegression
+from doubleml import DoubleMLLPQ, DoubleMLData
+from doubleml.datasets import make_iivm_data
+from doubleml.utils import dummy_classifier
+from ._utils import draw_smpls
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def normalize_ipw(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_lpq_fixture(dml_procedure, n_rep, normalize_ipw):
+    ext_predictions = {"d": {}}
+    np.random.seed(3141)
+    data = make_iivm_data(theta=0.5, n_obs=2000, dim_x=10, alpha_x=1.0, return_type="DataFrame")
+
+    dml_data = DoubleMLData(data, "y", "d", z_cols="z")
+    all_smpls = draw_smpls(len(dml_data.y), 5, n_rep=n_rep, groups=dml_data.d)
+
+    kwargs = {
+        "obj_dml_data": dml_data,
+        "score": "LPQ",
+        "n_rep": n_rep,
+        "dml_procedure": dml_procedure,
+        "normalize_ipw": normalize_ipw,
+        "draw_sample_splitting": False,
+    }
+
+    ml_g = LogisticRegression()
+    ml_m = LogisticRegression()
+
+    dml_lpq = DoubleMLLPQ(ml_g=ml_g, ml_m=ml_m, **kwargs)
+    dml_lpq.set_sample_splitting(all_smpls)
+
+    np.random.seed(3141)
+    dml_lpq.fit(store_predictions=True)
+
+    ext_predictions["d"]["ml_m_z"] = dml_lpq.predictions["ml_m_z"][:, :, 0]
+    ext_predictions["d"]["ml_m_d_z0"] = dml_lpq.predictions["ml_m_d_z0"][:, :, 0]
+    ext_predictions["d"]["ml_m_d_z1"] = dml_lpq.predictions["ml_m_d_z1"][:, :, 0]
+    ext_predictions["d"]["ml_g_du_z0"] = dml_lpq.predictions["ml_g_du_z0"][:, :, 0]
+    ext_predictions["d"]["ml_g_du_z1"] = dml_lpq.predictions["ml_g_du_z1"][:, :, 0]
+
+    dml_lpq_ext = DoubleMLLPQ(ml_g=dummy_classifier(), ml_m=dummy_classifier(), **kwargs)
+    dml_lpq_ext.set_sample_splitting(all_smpls)
+
+    np.random.seed(3141)
+    dml_lpq_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_lpq.coef[0], "coef_ext": dml_lpq_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_lpq_coef(doubleml_lpq_fixture):
+    assert math.isclose(doubleml_lpq_fixture["coef_normal"], doubleml_lpq_fixture["coef_ext"], rel_tol=1e-9, abs_tol=1e-4)
diff --git a/doubleml/tests/test_nonlinear_score_mixin.py b/doubleml/tests/test_nonlinear_score_mixin.py
index 2aca4f7a..ede9a0a1 100644
--- a/doubleml/tests/test_nonlinear_score_mixin.py
+++ b/doubleml/tests/test_nonlinear_score_mixin.py
@@ -100,7 +100,7 @@ def _check_score(self, score):
     def _check_data(self, obj_dml_data):
         pass
 
-    def _nuisance_est(self, smpls, n_jobs_cv, return_models=False):
+    def _nuisance_est(self, smpls, n_jobs_cv, external_predictions, return_models=False):
         x, y = check_X_y(self._dml_data.x, self._dml_data.y,
                          force_all_finite=False)
         x, d = check_X_y(x, self._dml_data.d,
diff --git a/doubleml/tests/test_pliv_external_predictions.py b/doubleml/tests/test_pliv_external_predictions.py
new file mode 100644
index 00000000..5a5eb097
--- /dev/null
+++ b/doubleml/tests/test_pliv_external_predictions.py
@@ -0,0 +1,100 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression
+from doubleml import DoubleMLPLIV, DoubleMLData
+from doubleml.datasets import make_pliv_CHS2015
+from doubleml.utils import dummy_regressor
+
+
+@pytest.fixture(scope="module", params=["partialling out", "IV-type"])
+def score(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def dim_z(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def adapted_doubleml_fixture(score, dml_procedure, n_rep, dim_z):
+    # IV-type score only allows dim_z = 1, so skip testcases with dim_z > 1 for IV-type score
+    if dim_z > 1 and score == "IV-type":
+        pytest.skip("IV-type score only allows dim_z = 1")
+        res_dict = None
+    else:
+        ext_predictions = {"d": {}}
+
+        data = make_pliv_CHS2015(
+            n_obs=500, dim_x=20, alpha=0.5, dim_z=dim_z, return_type="DataFrame"
+        )
+
+        np.random.seed(3141)
+
+        z_cols = [f"Z{i}" for i in range(1, dim_z + 1)]
+        dml_data = DoubleMLData(data, "y", "d", z_cols=z_cols)
+
+        kwargs = {
+            "obj_dml_data": dml_data,
+            "score": score,
+            "n_rep": n_rep,
+            "dml_procedure": dml_procedure,
+        }
+
+        if score == "IV-type":
+            kwargs["ml_g"] = LinearRegression()
+
+        dml_pliv = DoubleMLPLIV(
+            ml_m=LinearRegression(),
+            ml_l=LinearRegression(),
+            ml_r=LinearRegression(),
+            **kwargs,
+        )
+        np.random.seed(3141)
+
+        dml_pliv.fit(store_predictions=True)
+
+        ext_predictions["d"]["ml_l"] = dml_pliv.predictions["ml_l"][:, :, 0]
+        ext_predictions["d"]["ml_r"] = dml_pliv.predictions["ml_r"][:, :, 0]
+
+        if dim_z == 1:
+            ext_predictions["d"]["ml_m"] = dml_pliv.predictions["ml_m"][:, :, 0]
+            if score == "IV-type":
+                kwargs["ml_g"] = dummy_regressor()
+                ext_predictions["d"]["ml_g"] = dml_pliv.predictions["ml_g"][:, :, 0]
+        else:
+            for instr in range(dim_z):
+                ml_m_key = "ml_m_" + "Z" + str(instr + 1)
+                ext_predictions["d"][ml_m_key] = dml_pliv.predictions[ml_m_key][:, :, 0]
+
+        dml_pliv_ext = DoubleMLPLIV(
+            ml_m=dummy_regressor(), ml_l=dummy_regressor(), ml_r=dummy_regressor(), **kwargs
+        )
+
+        np.random.seed(3141)
+        dml_pliv_ext.fit(external_predictions=ext_predictions)
+
+        res_dict = {"coef_normal": dml_pliv.coef[0], "coef_ext": dml_pliv_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_adapted_doubleml_coef(adapted_doubleml_fixture):
+    assert math.isclose(
+        adapted_doubleml_fixture["coef_normal"],
+        adapted_doubleml_fixture["coef_ext"],
+        rel_tol=1e-9,
+        abs_tol=1e-4,
+    )
diff --git a/doubleml/tests/test_plr.py b/doubleml/tests/test_plr.py
index a27ee4a2..c3df97c4 100644
--- a/doubleml/tests/test_plr.py
+++ b/doubleml/tests/test_plr.py
@@ -77,10 +77,43 @@ def dml_plr_fixture(generate_data1, learner, score, dml_procedure):
     res_manual = fit_plr(y, x, d, clone(learner), clone(learner), clone(learner),
                          all_smpls, dml_procedure, score)
 
+    np.random.seed(3141)
+    # test with external nuisance predictions
+    if score == 'partialling out':
+        dml_plr_obj_ext = dml.DoubleMLPLR(obj_dml_data,
+                                      ml_l, ml_m,
+                                      n_folds,
+                                      score=score,
+                                      dml_procedure=dml_procedure)
+    else:
+        assert score == 'IV-type'
+        dml_plr_obj_ext = dml.DoubleMLPLR(obj_dml_data,
+                                      ml_l, ml_m, ml_g,
+                                      n_folds,
+                                      score=score,
+                                      dml_procedure=dml_procedure)
+
+    # synchronize the sample splitting
+    dml_plr_obj_ext.set_sample_splitting(all_smpls=all_smpls)
+
+    if score == 'partialling out':
+        prediction_dict = {'d': {'ml_l': dml_plr_obj.predictions['ml_l'].reshape(-1, 1),
+                                 'ml_m': dml_plr_obj.predictions['ml_m'].reshape(-1, 1)}}
+    else:
+        assert score == 'IV-type'
+        prediction_dict = {'d': {'ml_l': dml_plr_obj.predictions['ml_l'].reshape(-1, 1),
+                                 'ml_m': dml_plr_obj.predictions['ml_m'].reshape(-1, 1),
+                                 'ml_g': dml_plr_obj.predictions['ml_g'].reshape(-1, 1)}}
+
+    dml_plr_obj_ext.fit(external_predictions=prediction_dict)
+
+
     res_dict = {'coef': dml_plr_obj.coef,
                 'coef_manual': res_manual['theta'],
+                'coef_ext': dml_plr_obj_ext.coef,
                 'se': dml_plr_obj.se,
                 'se_manual': res_manual['se'],
+                'se_ext': dml_plr_obj_ext.se,
                 'boot_methods': boot_methods}
 
     for bootstrap in boot_methods:
@@ -91,10 +124,14 @@ def dml_plr_fixture(generate_data1, learner, score, dml_procedure):
 
         np.random.seed(3141)
         dml_plr_obj.bootstrap(method=bootstrap, n_rep_boot=n_rep_boot)
+        np.random.seed(3141)
+        dml_plr_obj_ext.bootstrap(method=bootstrap, n_rep_boot=n_rep_boot)
         res_dict['boot_coef' + bootstrap] = dml_plr_obj.boot_coef
         res_dict['boot_t_stat' + bootstrap] = dml_plr_obj.boot_t_stat
         res_dict['boot_coef' + bootstrap + '_manual'] = boot_theta
         res_dict['boot_t_stat' + bootstrap + '_manual'] = boot_t_stat
+        res_dict['boot_coef' + bootstrap + '_ext'] = dml_plr_obj_ext.boot_coef
+        res_dict['boot_t_stat' + bootstrap + '_ext'] = dml_plr_obj_ext.boot_t_stat
 
     # sensitivity tests
     res_dict['sensitivity_elements'] = dml_plr_obj.sensitivity_elements
@@ -114,6 +151,9 @@ def test_dml_plr_coef(dml_plr_fixture):
     assert math.isclose(dml_plr_fixture['coef'],
                         dml_plr_fixture['coef_manual'],
                         rel_tol=1e-9, abs_tol=1e-4)
+    assert math.isclose(dml_plr_fixture['coef'],
+                        dml_plr_fixture['coef_ext'],
+                        rel_tol=1e-9, abs_tol=1e-4)
 
 
 @pytest.mark.ci
@@ -121,6 +161,9 @@ def test_dml_plr_se(dml_plr_fixture):
     assert math.isclose(dml_plr_fixture['se'],
                         dml_plr_fixture['se_manual'],
                         rel_tol=1e-9, abs_tol=1e-4)
+    assert math.isclose(dml_plr_fixture['se'],
+                        dml_plr_fixture['se_ext'],
+                        rel_tol=1e-9, abs_tol=1e-4)
 
 
 @pytest.mark.ci
@@ -129,9 +172,15 @@ def test_dml_plr_boot(dml_plr_fixture):
         assert np.allclose(dml_plr_fixture['boot_coef' + bootstrap],
                            dml_plr_fixture['boot_coef' + bootstrap + '_manual'],
                            rtol=1e-9, atol=1e-4)
+        assert np.allclose(dml_plr_fixture['boot_coef' + bootstrap],
+                           dml_plr_fixture['boot_coef' + bootstrap + '_ext'],
+                           rtol=1e-9, atol=1e-4)
         assert np.allclose(dml_plr_fixture['boot_t_stat' + bootstrap],
                            dml_plr_fixture['boot_t_stat' + bootstrap + '_manual'],
                            rtol=1e-9, atol=1e-4)
+        assert np.allclose(dml_plr_fixture['boot_t_stat' + bootstrap],
+                           dml_plr_fixture['boot_t_stat' + bootstrap + '_ext'],
+                           rtol=1e-9, atol=1e-4)
 
 
 @pytest.mark.ci
diff --git a/doubleml/tests/test_plr_external_predictions.py b/doubleml/tests/test_plr_external_predictions.py
new file mode 100644
index 00000000..9d776710
--- /dev/null
+++ b/doubleml/tests/test_plr_external_predictions.py
@@ -0,0 +1,92 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LinearRegression
+from doubleml import DoubleMLPLR, DoubleMLData
+from doubleml.datasets import make_plr_CCDDHNR2018
+from doubleml.utils import dummy_regressor
+
+
+@pytest.fixture(scope="module", params=["IV-type", "partialling out"])
+def plr_score(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_m_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_l_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_g_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_plr_fixture(plr_score, dml_procedure, n_rep, set_ml_m_ext, set_ml_l_ext, set_ml_g_ext):
+    ext_predictions = {"d": {}}
+
+    x, y, d = make_plr_CCDDHNR2018(n_obs=500, dim_x=20, alpha=0.5, return_type="np.array")
+
+    np.random.seed(3141)
+
+    dml_data = DoubleMLData.from_arrays(x=x, y=y, d=d)
+
+    kwargs = {"obj_dml_data": dml_data, "score": plr_score, "n_rep": n_rep, "dml_procedure": dml_procedure}
+
+    if plr_score == "IV-type":
+        kwargs["ml_g"] = LinearRegression()
+
+    dml_plr = DoubleMLPLR(ml_m=LinearRegression(), ml_l=LinearRegression(), **kwargs)
+    np.random.seed(3141)
+
+    dml_plr.fit(store_predictions=True)
+
+    if set_ml_m_ext:
+        ext_predictions["d"]["ml_m"] = dml_plr.predictions["ml_m"][:, :, 0]
+        ml_m = dummy_regressor()
+    else:
+        ml_m = LinearRegression()
+
+    if set_ml_l_ext:
+        ext_predictions["d"]["ml_l"] = dml_plr.predictions["ml_l"][:, :, 0]
+        ml_l = dummy_regressor()
+    else:
+        ml_l = LinearRegression()
+
+    if plr_score == "IV-type" and set_ml_g_ext:
+        ext_predictions["d"]["ml_g"] = dml_plr.predictions["ml_g"][:, :, 0]
+        kwargs["ml_g"] = dummy_regressor()
+    elif plr_score == "IV-type" and not set_ml_g_ext:
+        kwargs["ml_g"] = LinearRegression()
+    else:
+        pass
+
+    dml_plr_ext = DoubleMLPLR(ml_m=ml_m, ml_l=ml_l, **kwargs)
+
+    np.random.seed(3141)
+    dml_plr_ext.fit(external_predictions=ext_predictions)
+
+    res_dict = {"coef_normal": dml_plr.coef[0], "coef_ext": dml_plr_ext.coef[0]}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_plr_coef(doubleml_plr_fixture):
+    assert math.isclose(doubleml_plr_fixture["coef_normal"], doubleml_plr_fixture["coef_ext"], rel_tol=1e-9, abs_tol=1e-4)
diff --git a/doubleml/tests/test_plr_rep_cross.py b/doubleml/tests/test_plr_rep_cross.py
index f2a50e21..9bbc2616 100644
--- a/doubleml/tests/test_plr_rep_cross.py
+++ b/doubleml/tests/test_plr_rep_cross.py
@@ -75,12 +75,45 @@ def dml_plr_fixture(generate_data1, learner, score, dml_procedure, n_rep):
     res_manual = fit_plr(y, x, d, _clone(learner), _clone(learner), _clone(learner),
                          all_smpls, dml_procedure, score, n_rep)
 
+    np.random.seed(3141)
+    # test with external nuisance predictions
+    if score == 'partialling out':
+        dml_plr_obj_ext = dml.DoubleMLPLR(obj_dml_data,
+                                      ml_l, ml_m,
+                                      n_folds,
+                                      n_rep,
+                                      score=score,
+                                      dml_procedure=dml_procedure)
+    else:
+        assert score == 'IV-type'
+        dml_plr_obj_ext = dml.DoubleMLPLR(obj_dml_data,
+                                      ml_l, ml_m, ml_g,
+                                      n_folds,
+                                      n_rep,
+                                      score=score,
+                                      dml_procedure=dml_procedure)
+
+    # synchronize the sample splitting
+    dml_plr_obj_ext.set_sample_splitting(all_smpls=all_smpls)
+
+    if score == 'partialling out':
+        prediction_dict = {'d': {'ml_l': dml_plr_obj.predictions['ml_l'].reshape(-1, n_rep),
+                                 'ml_m': dml_plr_obj.predictions['ml_m'].reshape(-1, n_rep)}}
+    else:
+        assert score == 'IV-type'
+        prediction_dict = {'d': {'ml_l': dml_plr_obj.predictions['ml_l'].reshape(-1, n_rep),
+                                 'ml_m': dml_plr_obj.predictions['ml_m'].reshape(-1, n_rep),
+                                 'ml_g': dml_plr_obj.predictions['ml_g'].reshape(-1, n_rep)}}
+
+    dml_plr_obj_ext.fit(external_predictions=prediction_dict)
+
     res_dict = {'coef': dml_plr_obj.coef,
                 'coef_manual': res_manual['theta'],
+                'coef_ext': dml_plr_obj_ext.coef,
                 'se': dml_plr_obj.se,
                 'se_manual': res_manual['se'],
-                'boot_methods': boot_methods
-                }
+                'se_ext': dml_plr_obj_ext.se,
+                'boot_methods': boot_methods}
 
     for bootstrap in boot_methods:
         np.random.seed(3141)
diff --git a/doubleml/tests/test_pq_external_predictions.py b/doubleml/tests/test_pq_external_predictions.py
new file mode 100644
index 00000000..4f23800f
--- /dev/null
+++ b/doubleml/tests/test_pq_external_predictions.py
@@ -0,0 +1,101 @@
+import numpy as np
+import pytest
+import math
+from sklearn.linear_model import LogisticRegression
+from doubleml import DoubleMLPQ, DoubleMLData
+from doubleml.datasets import make_irm_data
+from doubleml.utils import dummy_classifier
+from ._utils import draw_smpls
+
+
+@pytest.fixture(scope="module", params=["dml1", "dml2"])
+def dml_procedure(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[1, 3])
+def n_rep(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def normalize_ipw(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_m_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module", params=[True, False])
+def set_ml_g_ext(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def doubleml_pq_fixture(dml_procedure, n_rep, normalize_ipw, set_ml_m_ext, set_ml_g_ext):
+    ext_predictions = {"d": {}}
+    np.random.seed(3141)
+    data = make_irm_data(theta=1, n_obs=500, dim_x=5, return_type="DataFrame")
+
+    dml_data = DoubleMLData(data, "y", "d")
+    all_smpls = draw_smpls(len(dml_data.y), 5, n_rep=n_rep, groups=None)
+
+    kwargs = {
+        "obj_dml_data": dml_data,
+        "score": "PQ",
+        "n_rep": n_rep,
+        "dml_procedure": dml_procedure,
+        "normalize_ipw": normalize_ipw,
+        "draw_sample_splitting": False,
+    }
+
+    ml_m = LogisticRegression(random_state=42)
+    ml_g = LogisticRegression(random_state=42)
+
+    dml_pq = DoubleMLPQ(ml_g=ml_g, ml_m=ml_m, **kwargs)
+    dml_pq.set_sample_splitting(all_smpls)
+    np.random.seed(3141)
+    dml_pq.fit(store_predictions=True)
+
+    if set_ml_m_ext:
+        ext_predictions["d"]["ml_m"] = dml_pq.predictions["ml_m"][:, :, 0]
+        ml_m = dummy_classifier()
+    else:
+        ml_m = LogisticRegression(random_state=42)
+
+    if set_ml_g_ext:
+        ext_predictions["d"]["ml_g"] = dml_pq.predictions["ml_g"][:, :, 0]
+        ml_g = dummy_classifier()
+    else:
+        ml_g = LogisticRegression(random_state=42)
+
+    dml_pq_ext = DoubleMLPQ(ml_g=ml_g, ml_m=ml_m, **kwargs)
+    dml_pq_ext.set_sample_splitting(all_smpls)
+
+    np.random.seed(3141)
+    dml_pq_ext.fit(external_predictions=ext_predictions)
+
+    if set_ml_m_ext and not set_ml_g_ext:
+        # adjust tolerance for the case that ml_m is set to external predictions
+        # because no preliminary results are available for ml_m, the model use the (external) final predictions for ml_m
+        tol_rel = 0.1
+        tol_abs = 0.1
+    else:
+        tol_rel = 1e-9
+        tol_abs = 1e-4
+
+    res_dict = {"coef_normal": dml_pq.coef[0], "coef_ext": dml_pq_ext.coef[0], "tol_rel": tol_rel, "tol_abs": tol_abs}
+
+    return res_dict
+
+
+@pytest.mark.ci
+def test_doubleml_pq_coef(doubleml_pq_fixture):
+    assert math.isclose(
+        doubleml_pq_fixture["coef_normal"],
+        doubleml_pq_fixture["coef_ext"],
+        rel_tol=doubleml_pq_fixture["tol_rel"],
+        abs_tol=doubleml_pq_fixture["tol_abs"],
+    )