diff --git a/pyFTS/benchmarks/ResidualAnalysis.py b/pyFTS/benchmarks/ResidualAnalysis.py
index ae55eb5..fbfa16a 100644
--- a/pyFTS/benchmarks/ResidualAnalysis.py
+++ b/pyFTS/benchmarks/ResidualAnalysis.py
@@ -97,15 +97,17 @@ def plotResiduals(targets, models, tam=[8, 8], save=False, file=None):
 
 
 def plot_residuals(targets, models, tam=[8, 8], save=False, file=None):
-    fig, axes = plt.subplots(nrows=len(models), ncols=3, figsize=tam)
+    import scipy as sp
+    
+    fig, axes = plt.subplots(nrows=len(models), ncols=4, figsize=tam)
 
     for c, mfts in enumerate(models, start=0):
         if len(models) > 1:
             ax = axes[c]
         else:
             ax = axes
-        forecasts = mfts.forecast(targets)
-        res = residuals(targets, forecasts, mfts.order)
+        forecasts = mfts.predict(targets)
+        res = residuals(targets, forecasts, mfts.order+1)
         mu = np.mean(res)
         sig = np.std(res)
 
@@ -114,40 +116,46 @@ def plot_residuals(targets, models, tam=[8, 8], save=False, file=None):
         ax[0].set_xlabel(' ')
         ax[0].plot(res)
 
-        if c == 0: ax[1].set_title("Residuals Autocorrelation", size='large')
+        if c == 0: ax[1].set_title("Autocorrelation", size='large')
         ax[1].set_ylabel('ACS')
         ax[1].set_xlabel('Lag')
         ax[1].acorr(res)
 
-        if c == 0: ax[2].set_title("Residuals Histogram", size='large')
+        if c == 0: ax[2].set_title("Histogram", size='large')
         ax[2].set_ylabel('Freq')
         ax[2].set_xlabel('Bins')
         ax[2].hist(res)
+        
+        if c == 0: ax[3].set_title("QQ Plot", size='large')
+        
+        _, (__, ___, r) = sp.stats.probplot(res, plot=ax[3], fit=True)
 
     plt.tight_layout()
 
     Util.show_and_save_image(fig, file, save)
 
 
-def single_plot_residuals(targets, forecasts, order, tam=[8, 8], save=False, file=None):
-    fig, ax = plt.subplots(nrows=1, ncols=3, figsize=tam)
+def single_plot_residuals(targets, forecasts, order, tam=[10, 7], save=False, file=None):
+    import scipy as sp
+    
+    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=tam)
 
     res = residuals(targets, forecasts, order)
 
-    ax[0].set_title("Residuals", size='large')
-    ax[0].set_ylabel("Model", size='large')
-    ax[0].set_xlabel(' ')
-    ax[0].plot(res)
+    ax[0][0].set_title("Residuals", size='large')
+    ax[0][0].plot(res)
 
-    ax[1].set_title("Residuals Autocorrelation", size='large')
-    ax[1].set_ylabel('ACS')
-    ax[1].set_xlabel('Lag')
-    ax[1].acorr(res)
+    ax[0][1].set_title("Autocorrelation", size='large')
+    ax[0][1].set_ylabel('ACF')
+    ax[0][1].set_xlabel('Lag')
+    ax[0][1].acorr(res)
 
-    ax[2].set_title("Residuals Histogram", size='large')
-    ax[2].set_ylabel('Freq')
-    ax[2].set_xlabel('Bins')
-    ax[2].hist(res)
+    ax[1][0].set_title("Histogram", size='large')
+    ax[1][0].set_ylabel('Freq')
+    ax[1][0].set_xlabel('Bins')
+    ax[1][0].hist(res)
+    
+    _, (__, ___, r) = sp.stats.probplot(res, plot=ax[1][1], fit=True)
 
     plt.tight_layout()