diff --git a/tsfeatures/_modidx.py b/tsfeatures/_modidx.py
index c58ed29..f657af3 100644
--- a/tsfeatures/_modidx.py
+++ b/tsfeatures/_modidx.py
@@ -5,7 +5,8 @@
                 'doc_host': 'https://Nixtla.github.io',
                 'git_url': 'https://github.com/Nixtla/tsfeatures',
                 'lib_path': 'tsfeatures'},
-  'syms': { 'tsfeatures.m4_data': { 'tsfeatures.m4_data.m4_parser': ('m4_data.html#m4_parser', 'tsfeatures/m4_data.py'),
+  'syms': { 'tsfeatures.compare_with_r': {},
+            'tsfeatures.m4_data': { 'tsfeatures.m4_data.m4_parser': ('m4_data.html#m4_parser', 'tsfeatures/m4_data.py'),
                                     'tsfeatures.m4_data.maybe_download': ('m4_data.html#maybe_download', 'tsfeatures/m4_data.py'),
                                     'tsfeatures.m4_data.prepare_m4_data': ('m4_data.html#prepare_m4_data', 'tsfeatures/m4_data.py')},
             'tsfeatures.tsfeatures': { 'tsfeatures.tsfeatures._get_feats': ('tsfeatures.html#_get_feats', 'tsfeatures/tsfeatures.py'),
@@ -34,6 +35,7 @@
                                        'tsfeatures.tsfeatures.tsfeatures': ('tsfeatures.html#tsfeatures', 'tsfeatures/tsfeatures.py'),
                                        'tsfeatures.tsfeatures.unitroot_kpss': ('tsfeatures.html#unitroot_kpss', 'tsfeatures/tsfeatures.py'),
                                        'tsfeatures.tsfeatures.unitroot_pp': ('tsfeatures.html#unitroot_pp', 'tsfeatures/tsfeatures.py')},
+            'tsfeatures.tsfeatures_r': {},
             'tsfeatures.utils': { 'tsfeatures.utils.embed': ('utils.html#embed', 'tsfeatures/utils.py'),
                                   'tsfeatures.utils.hurst_exponent': ('utils.html#hurst_exponent', 'tsfeatures/utils.py'),
                                   'tsfeatures.utils.lambda_coef_var': ('utils.html#lambda_coef_var', 'tsfeatures/utils.py'),
diff --git a/tsfeatures/tsfeatures.py b/tsfeatures/tsfeatures.py
index 8c1d3cb..e905849 100644
--- a/tsfeatures/tsfeatures.py
+++ b/tsfeatures/tsfeatures.py
@@ -1,10 +1,10 @@
 # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/tsfeatures.ipynb.
 
 # %% auto 0
-__all__ = ['FREQS', 'acf_features', 'arch_stat', 'count_entropy', 'crossing_points', 'entropy', 'flat_spots', 'frequency',
-           'guerrero', 'heterogeneity', 'holt_parameters', 'hurst', 'hw_parameters', 'intervals', 'lumpiness',
-           'nonlinearity', 'pacf_features', 'series_length', 'sparsity', 'stability', 'stl_features', 'unitroot_kpss',
-           'unitroot_pp', 'tsfeatures']
+__all__ = ['acf_features', 'arch_stat', 'count_entropy', 'crossing_points', 'entropy', 'flat_spots', 'frequency', 'guerrero',
+           'heterogeneity', 'holt_parameters', 'hurst', 'hw_parameters', 'intervals', 'lumpiness', 'nonlinearity',
+           'pacf_features', 'series_length', 'sparsity', 'stability', 'stl_features', 'unitroot_kpss', 'unitroot_pp',
+           'tsfeatures']
 
 # %% ../nbs/tsfeatures.ipynb 3
 import os
@@ -45,9 +45,6 @@
 from .utils import *
 
 # %% ../nbs/tsfeatures.ipynb 7
-FREQS = {"H": 24, "D": 1, "M": 12, "Q": 4, "W": 1, "Y": 1}
-
-# %% ../nbs/tsfeatures.ipynb 8
 def acf_features(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Calculates autocorrelation function features.
 
@@ -87,7 +84,12 @@ def acf_features(x: np.array, freq: int = 1) -> Dict[str, float]:
     else:
         acfdiff2x = [np.nan] * 2
     # first autocorrelation coefficient
-    acf_1 = acfx[1]
+
+    try:
+        acf_1 = acfx[1]
+    except:
+        acf_1 = np.nan
+
     # sum of squares of first 10 autocorrelation coefficients
     sum_of_sq_acf10 = np.sum((acfx[1:11]) ** 2) if size_x > 10 else np.nan
     # first autocorrelation ciefficient of differenced series
@@ -113,7 +115,7 @@ def acf_features(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return output
 
-# %% ../nbs/tsfeatures.ipynb 11
+# %% ../nbs/tsfeatures.ipynb 10
 def arch_stat(
     x: np.array, freq: int = 1, lags: int = 12, demean: bool = True
 ) -> Dict[str, float]:
@@ -134,7 +136,7 @@ def arch_stat(
     if len(x) <= lags + 1:
         return {"arch_lm": np.nan}
     if demean:
-        x -= np.mean(x)
+        x = x - np.mean(x)
 
     size_x = len(x)
     mat = embed(x**2, lags + 1)
@@ -148,7 +150,7 @@ def arch_stat(
 
     return {"arch_lm": r_squared}
 
-# %% ../nbs/tsfeatures.ipynb 14
+# %% ../nbs/tsfeatures.ipynb 13
 def count_entropy(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Count entropy.
 
@@ -169,7 +171,7 @@ def count_entropy(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"count_entropy": entropy}
 
-# %% ../nbs/tsfeatures.ipynb 15
+# %% ../nbs/tsfeatures.ipynb 14
 def crossing_points(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Crossing points.
 
@@ -194,7 +196,7 @@ def crossing_points(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"crossing_points": cross.sum()}
 
-# %% ../nbs/tsfeatures.ipynb 16
+# %% ../nbs/tsfeatures.ipynb 15
 def entropy(x: np.array, freq: int = 1, base: float = e) -> Dict[str, float]:
     """Calculates sample entropy.
 
@@ -218,7 +220,7 @@ def entropy(x: np.array, freq: int = 1, base: float = e) -> Dict[str, float]:
 
     return {"entropy": entropy}
 
-# %% ../nbs/tsfeatures.ipynb 17
+# %% ../nbs/tsfeatures.ipynb 16
 def flat_spots(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Flat spots.
 
@@ -240,8 +242,9 @@ def flat_spots(x: np.array, freq: int = 1) -> Dict[str, float]:
         return {"flat_spots": np.nan}
 
     rlex = np.array([sum(1 for i in g) for k, g in groupby(cutx)]).max()
+    return {"flat_spots": rlex}
 
-# %% ../nbs/tsfeatures.ipynb 18
+# %% ../nbs/tsfeatures.ipynb 17
 def frequency(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Frequency.
 
@@ -260,7 +263,7 @@ def frequency(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"frequency": freq}
 
-# %% ../nbs/tsfeatures.ipynb 19
+# %% ../nbs/tsfeatures.ipynb 18
 def guerrero(
     x: np.array, freq: int = 1, lower: int = -1, upper: int = 2
 ) -> Dict[str, float]:
@@ -295,7 +298,7 @@ def guerrero(
 
     return {"guerrero": min_}
 
-# %% ../nbs/tsfeatures.ipynb 20
+# %% ../nbs/tsfeatures.ipynb 19
 def heterogeneity(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Heterogeneity.
 
@@ -358,7 +361,7 @@ def heterogeneity(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return output
 
-# %% ../nbs/tsfeatures.ipynb 21
+# %% ../nbs/tsfeatures.ipynb 20
 def holt_parameters(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Fitted parameters of a Holt model.
 
@@ -386,7 +389,7 @@ def holt_parameters(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return params
 
-# %% ../nbs/tsfeatures.ipynb 22
+# %% ../nbs/tsfeatures.ipynb 21
 def hurst(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Hurst index.
 
@@ -409,7 +412,7 @@ def hurst(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"hurst": hurst_index}
 
-# %% ../nbs/tsfeatures.ipynb 23
+# %% ../nbs/tsfeatures.ipynb 22
 def hw_parameters(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Fitted parameters of a Holt-Winters model.
 
@@ -441,7 +444,7 @@ def hw_parameters(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return params
 
-# %% ../nbs/tsfeatures.ipynb 24
+# %% ../nbs/tsfeatures.ipynb 23
 def intervals(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Intervals with demand.
 
@@ -465,7 +468,7 @@ def intervals(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"intervals_mean": np.mean(y), "intervals_sd": np.std(y, ddof=1)}
 
-# %% ../nbs/tsfeatures.ipynb 25
+# %% ../nbs/tsfeatures.ipynb 24
 def lumpiness(x: np.array, freq: int = 1) -> Dict[str, float]:
     """lumpiness.
 
@@ -499,7 +502,7 @@ def lumpiness(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"lumpiness": lumpiness}
 
-# %% ../nbs/tsfeatures.ipynb 26
+# %% ../nbs/tsfeatures.ipynb 25
 def nonlinearity(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Nonlinearity.
 
@@ -524,7 +527,7 @@ def nonlinearity(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"nonlinearity": test}
 
-# %% ../nbs/tsfeatures.ipynb 27
+# %% ../nbs/tsfeatures.ipynb 26
 def pacf_features(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Calculates partial autocorrelation function features.
 
@@ -595,7 +598,7 @@ def pacf_features(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return output
 
-# %% ../nbs/tsfeatures.ipynb 28
+# %% ../nbs/tsfeatures.ipynb 27
 def series_length(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Series length.
 
@@ -614,7 +617,7 @@ def series_length(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"series_length": len(x)}
 
-# %% ../nbs/tsfeatures.ipynb 29
+# %% ../nbs/tsfeatures.ipynb 28
 def sparsity(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Sparsity.
 
@@ -633,7 +636,7 @@ def sparsity(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"sparsity": np.mean(x == 0)}
 
-# %% ../nbs/tsfeatures.ipynb 30
+# %% ../nbs/tsfeatures.ipynb 29
 def stability(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Stability.
 
@@ -667,7 +670,7 @@ def stability(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"stability": stability}
 
-# %% ../nbs/tsfeatures.ipynb 31
+# %% ../nbs/tsfeatures.ipynb 30
 def stl_features(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Calculates seasonal trend using loess decomposition.
 
@@ -785,8 +788,14 @@ def stl_features(x: np.array, freq: int = 1) -> Dict[str, float]:
     time_x = add_constant(poly_m)
     coefs = OLS(trend0, time_x).fit().params
 
-    linearity = coefs[1]
-    curvature = -coefs[2]
+    try:
+        linearity = coefs[1]
+    except:
+        linearity = np.nan
+    try:
+        curvature = -coefs[2]
+    except:
+        curvature = np.nan
     # ACF features
     acfremainder = acf_features(remainder, m)
     # Assemble features
@@ -808,7 +817,7 @@ def stl_features(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return output
 
-# %% ../nbs/tsfeatures.ipynb 32
+# %% ../nbs/tsfeatures.ipynb 31
 def unitroot_kpss(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Unit root kpss.
 
@@ -834,7 +843,7 @@ def unitroot_kpss(x: np.array, freq: int = 1) -> Dict[str, float]:
 
     return {"unitroot_kpss": test_kpss}
 
-# %% ../nbs/tsfeatures.ipynb 33
+# %% ../nbs/tsfeatures.ipynb 32
 def unitroot_pp(x: np.array, freq: int = 1) -> Dict[str, float]:
     """Unit root pp.
 
@@ -884,7 +893,6 @@ def _get_feats(
     ],
     dict_freqs=FREQS,
 ):
-    print("dict_freq")
     if freq is None:
         inf_freq = pd.infer_freq(ts["ds"])
         if inf_freq is None: