diff --git a/toarstats/__init__.py b/toarstats/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..ff6d5cb3e771bca6bf49d0cb1e6854a1282f7f62 100644
--- a/toarstats/__init__.py
+++ b/toarstats/__init__.py
@@ -0,0 +1,10 @@
+"""The toarstats package.
+
+This package contains the following two subpackages:
+metrics - calculate metrics on time series
+trends - calculate trends on time series
+"""
+
+from toarstats import metrics, trends
+
+__all__ = ["metrics", "trends"]
diff --git a/toarstats/trends/__init__.py b/toarstats/trends/__init__.py
index 47d89bc5a3a0996228936a541fbe0514a1033801..052d97791348c418f5a641c8dbd4e2be8d8949c6 100644
--- a/toarstats/trends/__init__.py
+++ b/toarstats/trends/__init__.py
@@ -1,4 +1,12 @@
"""The trends subpackage.
This subpackage contains trend analysis tools for time series.
+
+This subpackage exports the following function:
+calculate_trend - public interface for the trends subpackage
"""
+
+from toarstats.trends.interface import calculate_trend
+
+
+__all__ = ["calculate_trend"]
diff --git a/toarstats/trends/interface.py b/toarstats/trends/interface.py
new file mode 100644
index 0000000000000000000000000000000000000000..f4ce6637549e9d0c46357000a36fd6a3cafc19a8
--- /dev/null
+++ b/toarstats/trends/interface.py
@@ -0,0 +1,46 @@
+"""This module contains the public interface for the trends subpackage.
+
+This module contains the following function:
+calculate_trend - calculate the trend using the requested method
+"""
+
+from toarstats.trends.ols import ols
+from toarstats.trends.quant_reg import quant_reg
+from toarstats.trends.utils import deseasonalize
+
+
+def calculate_trend(method, data, formula="value ~ datetime", quantiles=None):
+ """Calculate the trend using the requested method.
+
+ This function is the public interface for the ``trends`` subpackage.
+ It takes all the user inputs and returns the result of the requested
+ trend analysis.
+
+ :param method: either ``"OLS"`` or ``"quant"``
+ :param data: data containing a list of date time values and
+ associated parameter values on which to calculate the
+ trend
+ :param formula: the formula specifying the model
+ :param quantiles: a single quantile or a list of quantiles to
+ calculate, these must be between 0 and 1; only
+ needed when ``method="quant"``
+
+ :raises ValueError: raised if the ``method`` parameter is not
+ recognized
+
+ :return: The result of the fit or a list of fit results if
+ ``method="quant"`` and multiple quantiles are given
+ """
+ deseasonalized_data = deseasonalize(data)
+ if method == "OLS":
+ fit = ols(deseasonalized_data, formula)
+ elif method == "quant":
+ try:
+ fit = [quant_reg(deseasonalized_data, formula, quantile)
+ for quantile in quantiles]
+ except TypeError:
+ fit = quant_reg(deseasonalized_data, formula, quantiles)
+ else:
+ raise ValueError(f"{method} is not recognized, must be 'OLS' or"
+ " 'quant'")
+ return fit
diff --git a/toarstats/trends/ols.py b/toarstats/trends/ols.py
new file mode 100644
index 0000000000000000000000000000000000000000..255741d5a974bd450a2b8be0bdbc5bca972981b7
--- /dev/null
+++ b/toarstats/trends/ols.py
@@ -0,0 +1,22 @@
+"""Ordinary least squares (OLS) linear regression calculation.
+
+This module contains the following function:
+ols - calculate the OLS linear regression
+"""
+
+import statsmodels.formula.api as smf
+
+
+def ols(data, formula):
+ """Calculate the OLS linear regression.
+
+ :param data: data containing a list of date time values and
+ associated parameter values on which to calculate the
+ trend
+ :param formula: the formula specifying the model
+
+ :return: The result of the fit
+ """
+ model = smf.ols(formula, data)
+ fit = model.fit()
+ return fit
diff --git a/toarstats/trends/quant_reg.py b/toarstats/trends/quant_reg.py
new file mode 100644
index 0000000000000000000000000000000000000000..033597761266c22a61ad45aeece6b1212e0d00d6
--- /dev/null
+++ b/toarstats/trends/quant_reg.py
@@ -0,0 +1,23 @@
+"""Quantile regression calculation.
+
+This module contains the following function:
+quant_reg - calculate the OLS linear regression
+"""
+
+import statsmodels.formula.api as smf
+
+
+def quant_reg(data, formula, quantile):
+ """Calculate the quantile regression.
+
+ :param data: data containing a list of date time values and
+ associated parameter values on which to calculate the
+ trend
+ :param formula: the formula specifying the model
+ :param quantile: a single quantile, must be between 0 and 1
+
+ :return: The result of the fit
+ """
+ model = smf.quantreg(formula, data)
+ fit = model.fit(q=quantile)
+ return fit
diff --git a/toarstats/trends/utils.py b/toarstats/trends/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..0a07989d0c6fadc1115b6c3c4d0407ecbe982d2d
--- /dev/null
+++ b/toarstats/trends/utils.py
@@ -0,0 +1,29 @@
+"""This module contains helper functions for the trends subpackage.
+
+This module contains the following function:
+deseasonalize - calculate the trend using the requested method
+"""
+
+import numpy as np
+import pandas as pd
+import statsmodels.formula.api as smf
+
+
+def deseasonalize(data):
+ """Deseasonalize the data.
+
+ :param data: data containing a list of date time values and
+ associated parameter values
+
+ :return: The deseasonalized data
+ """
+ seasonal_cycle = smf.ols(
+ "value~np.cos(month*2*np.pi/12)+np.sin(month*2*np.pi/12)"
+ "+np.cos(month*4*np.pi/12)+np.sin(month*4*np.pi/12)",
+ pd.DataFrame({"value": data["value"], "month": data.index.month})
+ ).fit(method="qr").predict(pd.DataFrame({"month": range(1, 13)}))
+ idx = data.index
+ return pd.DataFrame({
+ "value": data["value"].values - seasonal_cycle[idx.month-1].values,
+ "datetime": (idx.year-2000)*12 + idx.month
+ })