diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py
index a253c63d..299d5e65 100644
--- a/GPy/models/__init__.py
+++ b/GPy/models/__init__.py
@@ -16,3 +16,4 @@ from gradient_checker import GradientChecker
 from ss_gplvm import SSGPLVM
 from gp_coregionalized_regression import GPCoregionalizedRegression
 from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression
+from gp_heteroscedastic_regression import GPHeteroscedasticRegression
diff --git a/GPy/models/gp_heteroscedastic_regression.py b/GPy/models/gp_heteroscedastic_regression.py
new file mode 100644
index 00000000..5a05fc98
--- /dev/null
+++ b/GPy/models/gp_heteroscedastic_regression.py
@@ -0,0 +1,39 @@
+# Copyright (c) 2012 - 2014 the GPy Austhors (see AUTHORS.txt)
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from ..core import GP
+from .. import likelihoods
+from .. import kern
+from .. import util
+
+class GPHeteroscedasticRegression(GP):
+    """
+    Gaussian Process model for heteroscedastic regression
+
+    This is a thin wrapper around the models.GP class, with a set of sensible defaults
+
+    :param X: input observations
+    :param Y: observed values
+    :param kernel: a GPy kernel, defaults to rbf
+    """
+    def __init__(self, X, Y, kernel=None, Y_metadata=None):
+
+        Ny = Y.shape[0]
+
+        if Y_metadata is None:
+            Y_metadata = {'output_index':np.arange(Ny)[:,None]}
+        else:
+            assert Y_metadata['output_index'].shape[0] == Ny
+
+        if kernel is None:
+            kernel = kern.RBF(X.shape[1])
+
+        #Likelihood
+        likelihoods_list = [likelihoods.Gaussian(name="Gaussian_noise_%s" %j) for j in range(Ny)]
+        likelihood = likelihoods.MixedNoise(likelihoods_list=likelihoods_list)
+
+        super(GPHeteroscedasticRegression, self).__init__(X,Y,kernel,likelihood, Y_metadata=Y_metadata)
+
+    def plot(self,*args):
+        return NotImplementedError