diff --git a/peakshaving_analyzer/input.py b/peakshaving_analyzer/input.py
index c437805..f1e5d12 100644
--- a/peakshaving_analyzer/input.py
+++ b/peakshaving_analyzer/input.py
@@ -87,6 +87,30 @@ def timeseries_to_df(self):
 
         return df
 
+    def calculate_statistics(self, print: bool = False) -> dict[str, float]:
+        ts_df = self.timeseries_to_df()
+        stats = {}
+
+        stats["min_load_kw"] = ts_df["consumption_kw"].min()
+        stats["max_load_kw"] = ts_df["consumption_kw"].max()
+        stats["mean_load_kw"] = ts_df["consumption_kw"].mean()
+        stats["median_load_kw"] = ts_df["consumption_kw"].median()
+        stats["variance"] = ts_df["consumption_kw"].var()
+        stats["std"] = ts_df["consumption_kw"].std()
+        stats["total_consumption_kwh"] = ts_df["consumption_kw"].sum() * self.hours_per_timestep
+
+        if print:
+            for key, value in stats.items():
+                print(f"{key}: {value}")
+
+        return stats
+
+    def plot_analysis(self):
+        self.plot_load_box()
+        self.plot_load_histogram()
+        self.plot_load_duration_curve()
+        self.plot_seasonal_decompose()
+
     def plot_load_duration_curve(self):
         ts_df = self.timeseries_to_df()
 
diff --git a/pyproject.toml b/pyproject.toml
index 701a01d..f996d59 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -37,6 +37,7 @@ dependencies = [
     "fine >=2.4.1",
     "pyomo >=6.8.0",
     "pgeocode >=0.5.0",
+    "statsmodels >=0.14.5",
     "highspy",
     "gurobipy",
     "argcomplete",