diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py
index 104a4fcfeda1..ac2a80b1d29d 100644
--- a/dpnp/dpnp_iface_statistics.py
+++ b/dpnp/dpnp_iface_statistics.py
@@ -54,7 +54,7 @@
     to_supported_dtypes,
 )
 
-from .dpnp_utils import call_origin, get_usm_allocations
+from .dpnp_utils import get_usm_allocations
 from .dpnp_utils.dpnp_utils_reduction import dpnp_wrap_reduction_call
 from .dpnp_utils.dpnp_utils_statistics import dpnp_cov, dpnp_median
 
@@ -748,61 +748,175 @@ def cov(
 
     For full documentation refer to :obj:`numpy.cov`.
 
+    Parameters
+    ----------
+    m : {dpnp.ndarray, usm_ndarray}
+        A 1-D or 2-D array containing multiple variables and observations.
+        Each row of `m` represents a variable, and each column a single
+        observation of all those variables. Also see `rowvar` below.
+    y : {None, dpnp.ndarray, usm_ndarray}, optional
+        An additional set of variables and observations. `y` has the same form
+        as that of `m`.
+
+        Default: ``None``.
+    rowvar : bool, optional
+        If `rowvar` is ``True``, then each row represents a variable, with
+        observations in the columns. Otherwise, the relationship is transposed:
+        each column represents a variable, while the rows contain observations.
+
+        Default: ``True``.
+    bias : bool, optional
+        Default normalization is by ``(N - 1)``, where ``N`` is the number of
+        observations given (unbiased estimate). If `bias` is ``True``, then
+        normalization is by ``N``. These values can be overridden by using the
+        keyword `ddof`.
+
+        Default: ``False``.
+    ddof : {None, int}, optional
+        If not ``None`` the default value implied by `bias` is overridden. Note
+        that ``ddof=1`` will return the unbiased estimate, even if both
+        `fweights` and `aweights` are specified, and ``ddof=0`` will return the
+        simple average. See the notes for the details.
+
+        Default: ``None``.
+    fweights : {None, dpnp.ndarray, usm_ndarray}, optional
+        1-D array of integer frequency weights; the number of times each
+        observation vector should be repeated.
+        It is required that ``fweights >= 0``. However, the function will not
+        raise an error when ``fweights < 0`` for performance reasons.
+
+        Default: ``None``.
+    aweights : {None, dpnp.ndarray, usm_ndarray}, optional
+        1-D array of observation vector weights. These relative weights are
+        typically large for observations considered "important" and smaller for
+        observations considered less "important". If ``ddof=0`` the array of
+        weights can be used to assign probabilities to observation vectors.
+        It is required that ``aweights >= 0``. However, the function will not
+        error when ``aweights < 0`` for performance reasons.
+
+        Default: ``None``.
+    dtype : {None, str, dtype object}, optional
+        Data-type of the result. By default, the return data-type will have
+        at least floating point type based on the capabilities of the device on
+        which the input arrays reside.
+
+        Default: ``None``.
+
     Returns
     -------
     out : dpnp.ndarray
         The covariance matrix of the variables.
 
-    Limitations
-    -----------
-    Input array ``m`` is supported as :obj:`dpnp.ndarray`.
-    Dimension of input array ``m`` is limited by ``m.ndim <= 2``.
-    Size and shape of input arrays are supported to be equal.
-    Parameter `y` is supported only with default value ``None``.
-    Parameter `bias` is supported only with default value ``False``.
-    Parameter `ddof` is supported only with default value ``None``.
-    Parameter `fweights` is supported only with default value ``None``.
-    Parameter `aweights` is supported only with default value ``None``.
-    Otherwise the function will be executed sequentially on CPU.
-    Input array data types are limited by supported DPNP :ref:`Data types`.
-
     See Also
     --------
-    :obj:`dpnp.corrcoef` : Normalized covariance matrix
+    :obj:`dpnp.corrcoef` : Normalized covariance matrix.
+
+    Notes
+    -----
+    Assume that the observations are in the columns of the observation array `m`
+    and let ``f = fweights`` and ``a = aweights`` for brevity. The steps to
+    compute the weighted covariance are as follows::
+
+        >>> import dpnp as np
+        >>> m = np.arange(10, dtype=np.float32)
+        >>> f = np.arange(10) * 2
+        >>> a = np.arange(10) ** 2.0
+        >>> ddof = 1
+        >>> w = f * a
+        >>> v1 = np.sum(w)
+        >>> v2 = np.sum(w * a)
+        >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
+        >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
+
+    Note that when ``a == 1``, the normalization factor
+    ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
+    as it should.
 
     Examples
     --------
     >>> import dpnp as np
     >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
-    >>> x.shape
-    (2, 3)
-    >>> [i for i in x]
-    [0, 1, 2, 2, 1, 0]
-    >>> out = np.cov(x)
-    >>> out.shape
-    (2, 2)
-    >>> [i for i in out]
-    [1.0, -1.0, -1.0, 1.0]
+
+    Consider two variables, :math:`x_0` and  :math:`x_1`, which correlate
+    perfectly, but in opposite directions:
+
+    >>> x
+    array([[0, 1, 2],
+           [2, 1, 0]])
+
+    Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
+    matrix shows this clearly:
+
+    >>> np.cov(x)
+    array([[ 1., -1.],
+           [-1.,  1.]])
+
+    Note that element :math:`C_{0,1}`, which shows the correlation between
+    :math:`x_0` and :math:`x_1`, is negative.
+
+    Further, note how `x` and `y` are combined:
+
+    >>> x = np.array([-2.1, -1,  4.3])
+    >>> y = np.array([3,  1.1,  0.12])
+    >>> X = np.stack((x, y), axis=0)
+    >>> np.cov(X)
+    array([[11.71      , -4.286     ], # may vary
+           [-4.286     ,  2.14413333]])
+    >>> np.cov(x, y)
+    array([[11.71      , -4.286     ], # may vary
+           [-4.286     ,  2.14413333]])
+    >>> np.cov(x)
+    array(11.71)
 
     """
 
-    if not dpnp.is_supported_array_type(m):
-        pass
-    elif m.ndim > 2:
-        pass
-    elif bias:
-        pass
-    elif ddof is not None:
-        pass
-    elif fweights is not None:
-        pass
-    elif aweights is not None:
-        pass
+    arrays = [m]
+    if y is not None:
+        arrays.append(y)
+    dpnp.check_supported_arrays_type(*arrays)
+
+    if m.ndim > 2:
+        raise ValueError("m has more than 2 dimensions")
+
+    if y is not None:
+        if y.ndim > 2:
+            raise ValueError("y has more than 2 dimensions")
+
+    if ddof is not None:
+        if not isinstance(ddof, int):
+            raise ValueError("ddof must be integer")
     else:
-        return dpnp_cov(m, y=y, rowvar=rowvar, dtype=dtype)
+        ddof = 0 if bias else 1
+
+    def_float = dpnp.default_float_type(m.sycl_queue)
+    if dtype is None:
+        dtype = dpnp.result_type(*arrays, def_float)
+
+    if fweights is not None:
+        dpnp.check_supported_arrays_type(fweights)
+        if not dpnp.issubdtype(fweights.dtype, numpy.integer):
+            raise TypeError("fweights must be integer")
+
+        if fweights.ndim > 1:
+            raise ValueError("cannot handle multidimensional fweights")
+
+        fweights = dpnp.astype(fweights, dtype=def_float)
+
+    if aweights is not None:
+        dpnp.check_supported_arrays_type(aweights)
+        if aweights.ndim > 1:
+            raise ValueError("cannot handle multidimensional aweights")
+
+        aweights = dpnp.astype(aweights, dtype=def_float)
 
-    return call_origin(
-        numpy.cov, m, y, rowvar, bias, ddof, fweights, aweights, dtype=dtype
+    return dpnp_cov(
+        m,
+        y=y,
+        rowvar=rowvar,
+        ddof=ddof,
+        dtype=dtype,
+        fweights=fweights,
+        aweights=aweights,
     )
 
 
diff --git a/dpnp/dpnp_utils/dpnp_utils_statistics.py b/dpnp/dpnp_utils/dpnp_utils_statistics.py
index 56cf5b011515..108fda7286fc 100644
--- a/dpnp/dpnp_utils/dpnp_utils_statistics.py
+++ b/dpnp/dpnp_utils/dpnp_utils_statistics.py
@@ -32,7 +32,6 @@
 
 import dpnp
 from dpnp.dpnp_array import dpnp_array
-from dpnp.dpnp_utils import get_usm_allocations, map_dtype_to_device
 
 __all__ = ["dpnp_cov", "dpnp_median"]
 
@@ -119,66 +118,77 @@ def _flatten_array_along_axes(a, axes_to_flatten, overwrite_input):
     return a_flatten, overwrite_input
 
 
-def dpnp_cov(m, y=None, rowvar=True, dtype=None):
+def dpnp_cov(
+    m, y=None, rowvar=True, ddof=1, dtype=None, fweights=None, aweights=None
+):
     """
-    dpnp_cov(m, y=None, rowvar=True, dtype=None)
-
     Estimate a covariance matrix based on passed data.
-    No support for given weights is provided now.
 
-    The implementation is done through existing dpnp and dpctl methods
-    instead of separate function call of dpnp backend.
+    The implementation is done through existing dpnp functions.
 
     """
 
-    def _get_2dmin_array(x, dtype):
-        """
-        Transform an input array to a form required for building a covariance matrix.
-
-        If applicable, it reshapes the input array to have 2 dimensions or greater.
-        If applicable, it transposes the input array when 'rowvar' is False.
-        It casts to another dtype, if the input array differs from requested one.
-
-        """
-        if x.ndim == 0:
-            x = x.reshape((1, 1))
-        elif x.ndim == 1:
-            x = x[dpnp.newaxis, :]
-
-        if not rowvar and x.ndim != 1:
-            x = x.T
-
-        if x.dtype != dtype:
-            x = dpnp.astype(x, dtype)
-        return x
+    # need to create a copy of input, since it will be modified in-place
+    x = dpnp.array(m, ndmin=2, dtype=dtype)
+    if not rowvar and m.ndim != 1:
+        x = x.T
 
-    # input arrays must follow CFD paradigm
-    _, queue = get_usm_allocations((m,) if y is None else (m, y))
+    if x.shape[0] == 0:
+        return dpnp.empty_like(
+            x, shape=(0, 0), dtype=dpnp.default_float_type(m.sycl_queue)
+        )
 
-    # calculate a type of result array if not passed explicitly
-    if dtype is None:
-        dtypes = [m.dtype, dpnp.default_float_type(sycl_queue=queue)]
-        if y is not None:
-            dtypes.append(y.dtype)
-        dtype = dpnp.result_type(*dtypes)
-        # TODO: remove when dpctl.result_type() is returned dtype based on fp64
-        dtype = map_dtype_to_device(dtype, queue.sycl_device)
-
-    X = _get_2dmin_array(m, dtype)
     if y is not None:
-        y = _get_2dmin_array(y, dtype)
-
-        X = dpnp.concatenate((X, y), axis=0)
-
-    avg = X.mean(axis=1)
+        y_ndim = y.ndim
+        y = dpnp.array(y, copy=None, ndmin=2, dtype=dtype)
+        if not rowvar and y_ndim != 1:
+            y = y.T
+        x = dpnp.concatenate((x, y), axis=0)
+
+    # get the product of frequencies and weights
+    w = None
+    if fweights is not None:
+        if fweights.shape[0] != x.shape[1]:
+            raise ValueError("incompatible numbers of samples and fweights")
+
+        w = fweights
+
+    if aweights is not None:
+        if aweights.shape[0] != x.shape[1]:
+            raise ValueError("incompatible numbers of samples and aweights")
+
+        if w is None:
+            w = aweights
+        else:
+            w *= aweights
+
+    avg, w_sum = dpnp.average(x, axis=1, weights=w, returned=True)
+    w_sum = w_sum[0]
+
+    # determine the normalization
+    if w is None:
+        fact = x.shape[1] - ddof
+    elif ddof == 0:
+        fact = w_sum
+    elif aweights is None:
+        fact = w_sum - ddof
+    else:
+        fact = w_sum - ddof * dpnp.sum(w * aweights) / w_sum
 
-    fact = X.shape[1] - 1
-    X -= avg[:, None]
+    if fact <= 0:
+        warnings.warn(
+            "Degrees of freedom <= 0 for slice", RuntimeWarning, stacklevel=2
+        )
+        fact = 0.0
 
-    c = dpnp.dot(X, X.T.conj())
-    c *= 1 / fact if fact != 0 else dpnp.nan
+    x -= avg[:, None]
+    if w is None:
+        x_t = x.T
+    else:
+        x_t = (x * w).T
 
-    return dpnp.squeeze(c)
+    c = dpnp.dot(x, x_t.conj()) / fact
+    return c.squeeze()
 
 
 def dpnp_median(
diff --git a/dpnp/tests/test_statistics.py b/dpnp/tests/test_statistics.py
index 0cc7b76c4498..93d7d7ea9277 100644
--- a/dpnp/tests/test_statistics.py
+++ b/dpnp/tests/test_statistics.py
@@ -5,6 +5,7 @@
 from numpy.testing import (
     assert_allclose,
     assert_array_equal,
+    assert_raises,
     assert_raises_regex,
 )
 
@@ -14,7 +15,9 @@
     assert_dtype_allclose,
     generate_random_numpy_array,
     get_all_dtypes,
+    get_complex_dtypes,
     get_float_complex_dtypes,
+    get_float_dtypes,
     has_support_aspect64,
     numpy_version,
 )
@@ -337,14 +340,236 @@ def test_correlate_unkown_method(self):
 
 class TestCov:
     @pytest.mark.parametrize(
-        "dtype", get_all_dtypes(no_bool=True, no_none=True, no_complex=True)
+        "dt", get_all_dtypes(no_none=True, no_complex=True)
     )
-    def test_false_rowvar_dtype(self, dtype):
-        a = numpy.array([[0, 2], [1, 1], [2, 0]], dtype=dtype)
+    def test_basic(self, dt):
+        a = numpy.array([[0, 2], [1, 1], [2, 0]], dtype=dt)
         ia = dpnp.array(a)
 
-        assert_allclose(dpnp.cov(ia.T), dpnp.cov(ia, rowvar=False))
-        assert_allclose(dpnp.cov(ia, rowvar=False), numpy.cov(a, rowvar=False))
+        expected = numpy.cov(a.T)
+        result = dpnp.cov(ia.T)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("dt", get_complex_dtypes())
+    def test_complex(self, dt):
+        a = numpy.array([[1, 2, 3], [1j, 2j, 3j]], dtype=dt)
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a)
+        result = dpnp.cov(ia)
+        assert_allclose(result, expected)
+
+        expected = numpy.cov(a, aweights=numpy.ones(3))
+        result = dpnp.cov(ia, aweights=dpnp.ones(3))
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize(
+        "dt", get_all_dtypes(no_none=True, no_complex=True)
+    )
+    @pytest.mark.parametrize("y_dt", get_complex_dtypes())
+    def test_y(self, dt, y_dt):
+        a = numpy.array([[1, 2, 3]], dtype=dt)
+        y = numpy.array([[1j, 2j, 3j]], dtype=y_dt)
+        ia, iy = dpnp.array(a), dpnp.array(y)
+
+        expected = numpy.cov(a, y)
+        result = dpnp.cov(ia, iy)
+        assert_allclose(result, expected)
+
+    @pytest.mark.filterwarnings("ignore::RuntimeWarning")
+    @pytest.mark.parametrize("sh", [None, (0, 2), (2, 0)])
+    def test_empty(self, sh):
+        a = numpy.array([]).reshape(sh)
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a)
+        result = dpnp.cov(ia)
+        assert_allclose(result, expected)
+
+    @pytest.mark.filterwarnings("ignore::RuntimeWarning")
+    def test_wrong_ddof(self):
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a.T, ddof=5)
+        result = dpnp.cov(ia.T, ddof=5)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("dt", get_float_dtypes())
+    @pytest.mark.parametrize("rowvar", [True, False])
+    def test_1D_rowvar(self, dt, rowvar):
+        a = numpy.array([0.3942, 0.5969, 0.7730, 0.9918, 0.7964], dtype=dt)
+        y = numpy.array([0.0780, 0.3107, 0.2111, 0.0334, 0.8501])
+        ia, iy = dpnp.array(a), dpnp.array(y)
+
+        expected = numpy.cov(a, rowvar=rowvar)
+        result = dpnp.cov(ia, rowvar=rowvar)
+        assert_allclose(result, expected)
+
+        expected = numpy.cov(a, y, rowvar=rowvar)
+        result = dpnp.cov(ia, iy, rowvar=rowvar)
+        assert_allclose(result, expected)
+
+    def test_1D_variance(self):
+        a = numpy.array([0.3942, 0.5969, 0.7730, 0.9918, 0.7964])
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a, ddof=1)
+        result = dpnp.cov(ia, ddof=1)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("freq_data", [[1, 4, 1], [1, 1, 1]])
+    def test_fweights(self, freq_data):
+        a = numpy.array([0.0, 1.0, 2.0], ndmin=2)
+        freq = numpy.array(freq_data)
+        ia, ifreq = dpnp.array(a), dpnp.array(freq_data)
+
+        expected = numpy.cov(a, fweights=freq)
+        result = dpnp.cov(ia, fweights=ifreq)
+        assert_allclose(result, expected)
+
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a.T, fweights=freq)
+        result = dpnp.cov(ia.T, fweights=ifreq)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("xp", [dpnp, numpy])
+    def test_float_fweights(self, xp):
+        a = xp.array([[0, 2], [1, 1], [2, 0]])
+        freq = xp.array([1, 4, 1]) + 0.5
+        assert_raises(TypeError, xp.cov, a, fweights=freq)
+
+    @pytest.mark.parametrize("xp", [dpnp, numpy])
+    @pytest.mark.parametrize("sh", [(2, 3), 2])
+    def test_fweights_wrong_shapes(self, xp, sh):
+        a = xp.array([[0, 2], [1, 1], [2, 0]])
+        freq = xp.ones(sh, dtype=xp.int_)
+        assert_raises((ValueError, RuntimeError), xp.cov, a.T, fweights=freq)
+
+    @pytest.mark.parametrize("freq", [numpy.array([1, 4, 1]), 2])
+    def test_fweights_wrong_type(self, freq):
+        a = dpnp.array([[0, 2], [1, 1], [2, 0]]).T
+        assert_raises(TypeError, dpnp.cov, a, fweights=freq)
+
+    @pytest.mark.parametrize("weights_data", [[1.0, 4.0, 1.0], [1.0, 1.0, 1.0]])
+    def test_aweights(self, weights_data):
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        weights = numpy.array(weights_data)
+        ia, iweights = dpnp.array(a), dpnp.array(weights_data)
+
+        expected = numpy.cov(a.T, aweights=weights)
+        result = dpnp.cov(ia.T, aweights=iweights)
+        assert_allclose(result, expected)
+
+        expected = numpy.cov(a.T, aweights=3.0 * weights)
+        result = dpnp.cov(ia.T, aweights=3.0 * iweights)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("xp", [dpnp, numpy])
+    @pytest.mark.parametrize("sh", [(2, 3), 2])
+    def test_aweights_wrong_shapes(self, xp, sh):
+        a = xp.array([[0, 2], [1, 1], [2, 0]])
+        weights = xp.ones(sh)
+        assert_raises((ValueError, RuntimeError), xp.cov, a.T, aweights=weights)
+
+    @pytest.mark.parametrize("weights", [numpy.array([1.0, 4.0, 1.0]), 2.0])
+    def test_aweights_wrong_type(self, weights):
+        a = dpnp.array([[0, 2], [1, 1], [2, 0]]).T
+        assert_raises(TypeError, dpnp.cov, a, aweights=weights)
+
+    def test_unit_fweights_and_aweights(self):
+        a = numpy.array([0.0, 1.0, 2.0], ndmin=2)
+        freq = numpy.array([1, 4, 1])
+        weights = numpy.ones(3)
+        ia, ifreq, iweights = (
+            dpnp.array(a),
+            dpnp.array(freq),
+            dpnp.array(weights),
+        )
+
+        # unit weights
+        expected = numpy.cov(a, fweights=freq, aweights=weights)
+        result = dpnp.cov(ia, fweights=ifreq, aweights=iweights)
+        assert_allclose(result, expected)
+
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        ia = dpnp.array(a)
+
+        # unit weights
+        expected = numpy.cov(a.T, fweights=freq, aweights=weights)
+        result = dpnp.cov(ia.T, fweights=ifreq, aweights=iweights)
+        assert_allclose(result, expected)
+
+        freq = numpy.ones(3, dtype=numpy.int_)
+        ifreq = dpnp.array(freq)
+
+        # unit frequencies and weights
+        expected = numpy.cov(a.T, fweights=freq, aweights=weights)
+        result = dpnp.cov(ia.T, fweights=ifreq, aweights=iweights)
+        assert_allclose(result, expected)
+
+        weights = numpy.array([1.0, 4.0, 1.0])
+        iweights = dpnp.array(weights)
+
+        # unit frequencies
+        expected = numpy.cov(a.T, fweights=freq, aweights=weights)
+        result = dpnp.cov(ia.T, fweights=ifreq, aweights=iweights)
+        assert_allclose(result, expected)
+
+        expected = numpy.cov(a.T, fweights=freq, aweights=3.0 * weights)
+        result = dpnp.cov(ia.T, fweights=ifreq, aweights=3.0 * iweights)
+        assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("dt", get_float_complex_dtypes())
+    def test_dtype(self, dt):
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a.T, dtype=dt)
+        result = dpnp.cov(ia.T, dtype=dt)
+        assert_allclose(result, expected)
+        assert result.dtype == dt
+
+    @pytest.mark.parametrize("dt", get_float_complex_dtypes())
+    @pytest.mark.parametrize("bias", [True, False])
+    def test_bias(self, dt, bias):
+        a = generate_random_numpy_array((3, 4), dtype=dt)
+        ia = dpnp.array(a)
+
+        expected = numpy.cov(a, bias=bias)
+        result = dpnp.cov(ia, bias=bias)
+        assert_dtype_allclose(result, expected)
+
+        # with rowvar
+        expected = numpy.cov(a, rowvar=False, bias=bias)
+        result = dpnp.cov(ia, rowvar=False, bias=bias)
+        assert_dtype_allclose(result, expected)
+
+        freq = numpy.array([1, 4, 1, 7])
+        ifreq = dpnp.array(freq)
+
+        # with frequency
+        expected = numpy.cov(a, bias=bias, fweights=freq)
+        result = dpnp.cov(ia, bias=bias, fweights=ifreq)
+        assert_dtype_allclose(result, expected)
+
+        weights = numpy.array([1.2, 3.7, 5.0, 1.1])
+        iweights = dpnp.array(weights)
+
+        # with weights
+        expected = numpy.cov(a, bias=bias, aweights=weights)
+        result = dpnp.cov(ia, bias=bias, aweights=iweights)
+        assert_dtype_allclose(result, expected)
+
+    def test_usm_ndarray(self):
+        a = numpy.array([[0, 2], [1, 1], [2, 0]])
+        ia = dpt.asarray(a)
+
+        expected = numpy.cov(a.T)
+        result = dpnp.cov(ia.T)
+        assert_allclose(result, expected)
 
     # numpy 2.2 properly transposes 2d array when rowvar=False
     @with_requires("numpy>=2.2")
diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py
index c1041760823d..b82df205746e 100644
--- a/dpnp/tests/test_sycl_queue.py
+++ b/dpnp/tests/test_sycl_queue.py
@@ -463,6 +463,7 @@ def test_meshgrid(device):
             "cos", [-dpnp.pi / 2, -dpnp.pi / 4, 0.0, dpnp.pi / 4, dpnp.pi / 2]
         ),
         pytest.param("cosh", [-5.0, -3.5, 0.0, 3.5, 5.0]),
+        pytest.param("cov", [[0, 1, 2], [2, 1, 0]]),
         pytest.param("count_nonzero", [3, 0, 2, -1.2]),
         pytest.param("cumprod", [[1, 2, 3], [4, 5, 6]]),
         pytest.param("cumsum", [[1, 2, 3], [4, 5, 6]]),
@@ -729,6 +730,7 @@ def test_reduce_hypot(device):
             [[0.7, 0.8, 0.9], [1.0, 1.1, 1.2]],
         ),
         pytest.param("correlate", [1, 2, 3], [4, 5, 6]),
+        pytest.param("cov", [-2.1, -1, 4.3], [3, 1.1, 0.12]),
         pytest.param("cross", [1.0, 2.0, 3.0], [4.0, 5.0, 6.0]),
         pytest.param("digitize", [0.2, 6.4, 3.0], [0.0, 1.0, 2.5, 4.0]),
         pytest.param(
diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py
index b3452e0ccf72..c0024409166c 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -621,6 +621,7 @@ def test_norm(usm_type, ord, axis):
             "cos", [-dp.pi / 2, -dp.pi / 4, 0.0, dp.pi / 4, dp.pi / 2]
         ),
         pytest.param("cosh", [-5.0, -3.5, 0.0, 3.5, 5.0]),
+        pytest.param("cov", [[0, 1, 2], [2, 1, 0]]),
         pytest.param("count_nonzero", [0, 1, 7, 0]),
         pytest.param("cumlogsumexp", [1.0, 2.0, 4.0, 7.0]),
         pytest.param("cumprod", [[1, 2, 3], [4, 5, 6]]),
@@ -737,6 +738,7 @@ def test_1in_1out(func, data, usm_type):
             [[0.7, 0.8, 0.9], [1.0, 1.1, 1.2]],
         ),
         pytest.param("correlate", [1, 2, 3], [0, 1, 0.5]),
+        pytest.param("cov", [-2.1, -1, 4.3], [3, 1.1, 0.12]),
         # dpnp.dot has 3 different implementations based on input arrays dtype
         # checking all of them
         pytest.param("dot", [3.0, 4.0, 5.0], [1.0, 2.0, 3.0]),
diff --git a/dpnp/tests/third_party/cupy/statistics_tests/test_correlation.py b/dpnp/tests/third_party/cupy/statistics_tests/test_correlation.py
index 66aeb374085c..604e545e0785 100644
--- a/dpnp/tests/third_party/cupy/statistics_tests/test_correlation.py
+++ b/dpnp/tests/third_party/cupy/statistics_tests/test_correlation.py
@@ -7,6 +7,11 @@
 from dpnp.tests.helper import has_support_aspect64, numpy_version
 from dpnp.tests.third_party.cupy import testing
 
+if numpy_version() >= "2.0.0":
+    from numpy._core._exceptions import _UFuncOutputCastingError
+else:
+    from numpy.core._exceptions import _UFuncOutputCastingError
+
 
 class TestCorrcoef(unittest.TestCase):
 
@@ -60,6 +65,26 @@ def generate_input(self, a_shape, y_shape, xp, dtype):
             y = testing.shaped_arange(y_shape, xp, dtype)
         return a, y
 
+    def call_cov(self, xp, a, y, rowvar, bias, ddof, fweights, aweights, dtype):
+        try:
+            return xp.cov(
+                a, y, rowvar, bias, ddof, fweights, aweights, dtype=dtype
+            )
+        except ValueError as e:
+            if (
+                xp is cupy
+                and "function 'subtract' does not support input types" in str(e)
+            ):
+                # numpy raises _UFuncOutputCastingError
+                raise _UFuncOutputCastingError(
+                    numpy.subtract,
+                    "same_kind",
+                    numpy.dtype("f8"),
+                    numpy.dtype(dtype),
+                    0,
+                )
+            raise
+
     @testing.for_all_dtypes()
     @testing.numpy_cupy_allclose(
         type_check=has_support_aspect64(), accept_error=True
@@ -82,10 +107,9 @@ def check(
             fweights = name.asarray(fweights)
         if aweights is not None:
             aweights = name.asarray(aweights)
-        # print(type(fweights))
-        # return xp.cov(a, y, rowvar, bias, ddof,
-        #               fweights, aweights, dtype=dtype)
-        return xp.cov(a, y, rowvar, bias, ddof, fweights, aweights)
+        return self.call_cov(
+            xp, a, y, rowvar, bias, ddof, fweights, aweights, dtype
+        )
 
     @testing.for_all_dtypes()
     @testing.numpy_cupy_allclose(accept_error=True)
@@ -103,8 +127,8 @@ def check_warns(
     ):
         with testing.assert_warns(RuntimeWarning):
             a, y = self.generate_input(a_shape, y_shape, xp, dtype)
-            return xp.cov(
-                a, y, rowvar, bias, ddof, fweights, aweights, dtype=dtype
+            return self.call_cov(
+                xp, a, y, rowvar, bias, ddof, fweights, aweights, dtype
             )
 
     @testing.for_all_dtypes()
@@ -126,8 +150,8 @@ def check_raises(
                     a, y, rowvar, bias, ddof, fweights, aweights, dtype=dtype
                 )
 
-    @pytest.mark.usefixtures("allow_fall_back_on_numpy")
     @pytest.mark.filterwarnings("ignore::RuntimeWarning")
+    # @testing.with_requires("numpy>=2.2")
     def test_cov(self):
         self.check((2, 3))
         self.check((2,), (2,))
@@ -144,12 +168,10 @@ def test_cov(self):
         self.check((1, 3), bias=True, aweights=(1.0, 4.0, 1.0))
         self.check((1, 3), fweights=(1, 4, 1), aweights=(1.0, 4.0, 1.0))
 
-    @pytest.mark.usefixtures("allow_fall_back_on_numpy")
     def test_cov_warns(self):
         self.check_warns((2, 3), ddof=3)
         self.check_warns((2, 3), ddof=4)
 
-    @pytest.mark.usefixtures("allow_fall_back_on_numpy")
     def test_cov_raises(self):
         self.check_raises((2, 3), ddof=1.2)
         self.check_raises((3, 4, 2))