diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 2f4b7db46b4..8dccec79378 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -9,6 +9,7 @@
 - Pass the `tune` argument from `sample` when using `advi+adapt_diag_grad` (see issue [#3965](https://github.com/pymc-devs/pymc3/issues/3965), fixed by [#3979](https://github.com/pymc-devs/pymc3/pull/3979)).
 - Add simple test case for new coords and dims feature in `pm.Model` (see [#3977](https://github.com/pymc-devs/pymc3/pull/3977)). 
 - Require ArviZ >= 0.9.0 (see [#3977](https://github.com/pymc-devs/pymc3/pull/3977)).
+- Fixed issue [#3962](https://github.com/pymc-devs/pymc3/issues/3962) by making change in the `_random()` method of `GaussianRandomWalk` class, refer to PR [#3985]. Further testing revealed a new issue which is being tracked [#4010](https://github.com/pymc-devs/pymc3/issues/4010) 
 
 _NB: The `docs/*` folder is still removed from the tarball due to an upload size limit on PyPi._
 
diff --git a/docs/release-notes/pymc3-3.0.md b/docs/release-notes/pymc3-3.0.md
index 8b137891791..e69de29bb2d 100644
--- a/docs/release-notes/pymc3-3.0.md
+++ b/docs/release-notes/pymc3-3.0.md
@@ -1 +0,0 @@
-
diff --git a/pymc3/distributions/timeseries.py b/pymc3/distributions/timeseries.py
index 56351a53c25..d5d836b4655 100644
--- a/pymc3/distributions/timeseries.py
+++ b/pymc3/distributions/timeseries.py
@@ -17,6 +17,7 @@
 from scipy import stats
 import theano.tensor as tt
 from theano import scan
+import numpy as np
 
 from pymc3.util import get_variable_name
 from .continuous import get_tau_sigma, Normal, Flat
@@ -303,14 +304,23 @@ def random(self, point=None, size=None):
         )
 
     def _random(self, sigma, mu, size, sample_shape):
-        """Implement a Gaussian random walk as a cumulative sum of normals."""
+        """Implement a Gaussian random walk as a cumulative sum of normals.
+        axis = len(size) - 1 denotes the axis along which cumulative sum would be calculated.
+        This might need to be corrected in future when issue #4010 is fixed.
+        Lines 318-322 ties the starting point of each instance of random walk to 0"
+        """
         if size[len(sample_shape)] == sample_shape:
             axis = len(sample_shape)
         else:
-            axis = 0
+            axis = len(size) - 1
         rv = stats.norm(mu, sigma)
         data = rv.rvs(size).cumsum(axis=axis)
-        data = data - data[0]  # TODO: this should be a draw from `init`, if available
+        data = np.array(data)
+        if len(data.shape)>1:
+            for i in range(data.shape[0]):
+                data[i] = data[i] - data[i][0]
+        else:
+            data = data - data[0]
         return data
 
     def _repr_latex_(self, name=None, dist=None):