diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index c4b4ce2..b136f64 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -1,23 +1,16 @@
 # Read the Docs configuration file
 # See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
 
-# Required
 version: 2
 
-# Set the OS, Python version, and other tools you might need
 build:
   os: ubuntu-24.04
   tools:
     python: "3.13"
 
-# Build documentation in the "docs/" directory with Sphinx
 sphinx:
    configuration: docs/conf.py
 
-# Optionally, but recommended,
-# declare the Python requirements required to build your documentation
-# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
-# python:
-#    install:
-#    - requirements: docs/requirements.txt
-        
\ No newline at end of file
+python:
+  install:
+    - requirements: docs/requirements.txt
diff --git a/README.md b/README.md
index cc40109..605fd82 100644
--- a/README.md
+++ b/README.md
@@ -11,30 +11,31 @@
 ## Installation
 
 ```bash
+pip install dynaris
+# or
 uv add dynaris
 ```
 
+## Documentation
+
+Full documentation is available at [dynaris.readthedocs.io](https://dynaris.readthedocs.io).
+
 ## Quickstart
 
 ```python
-from dynaris import LocalLevel, Seasonal, DLM
+from dynaris import LocalLevel, DLM
 from dynaris.datasets import load_nile
 
 # Load data
 y = load_nile()
 
-# Build a model by composing components
-model = LocalLevel(sigma_level=38.0, sigma_obs=123.0) + Seasonal(period=12)
-
-# Fit, smooth, forecast
-dlm = DLM(model)
+# Build a local-level model and fit
+dlm = DLM(LocalLevel(sigma_level=38.0, sigma_obs=123.0))
 dlm.fit(y).smooth()
-fc = dlm.forecast(steps=12)
 
-# Print summary
+# Forecast and plot
+fc = dlm.forecast(steps=10)
 print(dlm.summary())
-
-# Single-figure overview
 dlm.plot(kind="panel")
 ```
 
@@ -43,7 +44,7 @@ dlm.plot(kind="panel")
 Build models by combining components with `+`:
 
 ```python
-from dynaris import LocalLinearTrend, Seasonal, Cycle, Autoregressive, Regression
+from dynaris import LocalLinearTrend, Seasonal, Cycle
 
 model = (
     LocalLinearTrend(sigma_level=1.0, sigma_slope=0.1)
@@ -80,26 +81,14 @@ print(f"Log-likelihood: {result.log_likelihood:.2f}")
 
 ## Datasets
 
-```python
-from dynaris.datasets import load_nile, load_airline, load_lynx, load_sunspots, load_temperature, load_gdp
-
-y = load_airline()     # 144 monthly obs, 1949-1960
-y = load_lynx()        # 114 annual obs, 1821-1934 (~10-year cycle)
-y = load_sunspots()    # 288 annual obs, 1700-1987 (~11-year cycle)
-y = load_temperature() # 144 annual obs, 1880-2023 (warming trend)
-y = load_gdp()         # 319 quarterly obs, 1947-2026 (business cycle)
-```
-
-## Notation
-
-Dynaris follows the West & Harrison (1997) notation:
-
-| Symbol | Code | Meaning |
-|--------|------|---------|
-| **G** | `model.G` / `system_matrix` | System (evolution) matrix |
-| **F** | `model.F` / `observation_matrix` | Observation (regression) matrix |
-| **W** | `model.W` / `evolution_cov` | Evolution covariance |
-| **V** | `model.V` / `obs_cov` | Observational variance |
+| Dataset | Loader | N | Frequency | Domain |
+|---------|--------|---|-----------|--------|
+| Nile river flow | `load_nile()` | 100 | Annual | Hydrology |
+| Airline passengers | `load_airline()` | 144 | Monthly | Transportation |
+| Lynx population | `load_lynx()` | 114 | Annual | Ecology |
+| Sunspot numbers | `load_sunspots()` | 288 | Annual | Astronomy |
+| Global temperature | `load_temperature()` | 144 | Annual | Climate |
+| US GDP growth | `load_gdp()` | 319 | Quarterly | Economics |
 
 ## License
 
diff --git a/docs/_static/custom.css b/docs/_static/custom.css
new file mode 100644
index 0000000..23f074f
--- /dev/null
+++ b/docs/_static/custom.css
@@ -0,0 +1,3 @@
+.sidebar-logo img {
+    max-width: 80px;
+}
diff --git a/docs/_static/logo.png b/docs/_static/logo.png
new file mode 100644
index 0000000..9efcc20
Binary files /dev/null and b/docs/_static/logo.png differ
diff --git a/docs/api/components.rst b/docs/api/components.rst
index bfa239a..cecc63e 100644
--- a/docs/api/components.rst
+++ b/docs/api/components.rst
@@ -3,16 +3,25 @@ DLM Components
 
 Composable building blocks for Dynamic Linear Models. Each function returns
 a :class:`~dynaris.core.state_space.StateSpaceModel` that can be composed
-via the ``+`` operator.
+via the ``+`` operator. See :doc:`/user-guide/components` for usage guidance.
+
+Trend
+-----
 
 .. autofunction:: dynaris.dlm.components.LocalLevel
 
 .. autofunction:: dynaris.dlm.components.LocalLinearTrend
 
+Periodic
+--------
+
 .. autofunction:: dynaris.dlm.components.Seasonal
 
+.. autofunction:: dynaris.dlm.components.Cycle
+
+Other
+-----
+
 .. autofunction:: dynaris.dlm.components.Regression
 
 .. autofunction:: dynaris.dlm.components.Autoregressive
-
-.. autofunction:: dynaris.dlm.components.Cycle
diff --git a/docs/api/core.rst b/docs/api/core.rst
index ad5a522..5401284 100644
--- a/docs/api/core.rst
+++ b/docs/api/core.rst
@@ -1,12 +1,16 @@
 Core Types
 ==========
 
-Fundamental data structures: state-space model, Gaussian state, result
-containers, and filter/smoother protocols.
+Fundamental data structures used throughout dynaris. These are the building
+blocks that filters, smoothers, and the DLM class operate on.
 
 StateSpaceModel
 ---------------
 
+The central model representation. Holds the four system matrices (F, G, V, W)
+following West and Harrison (1997) notation. Returned by all component
+functions and composed via ``+``.
+
 .. autoclass:: dynaris.core.state_space.StateSpaceModel
    :members:
    :show-inheritance:
@@ -14,6 +18,9 @@ StateSpaceModel
 GaussianState
 -------------
 
+Represents a Gaussian belief about the state: a mean vector and covariance
+matrix. Used internally by the Kalman filter and smoother at each time step.
+
 .. autoclass:: dynaris.core.types.GaussianState
    :members:
    :show-inheritance:
@@ -21,20 +28,31 @@ GaussianState
 FilterResult
 ------------
 
+Container returned by the Kalman filter. Holds filtered state means,
+covariances, log-likelihood, and forecast errors for all time steps.
+
 .. autoclass:: dynaris.core.results.FilterResult
    :members:
    :show-inheritance:
+   :no-index:
 
 SmootherResult
 --------------
 
+Container returned by the RTS smoother. Holds smoothed state means and
+covariances for all time steps.
+
 .. autoclass:: dynaris.core.results.SmootherResult
    :members:
    :show-inheritance:
+   :no-index:
 
 Protocols
 ---------
 
+Interfaces that filter and smoother implementations must satisfy. Useful
+for type checking and extending dynaris with custom algorithms.
+
 .. autoclass:: dynaris.core.protocols.FilterProtocol
    :members:
 
diff --git a/docs/api/datasets.rst b/docs/api/datasets.rst
new file mode 100644
index 0000000..73fcd17
--- /dev/null
+++ b/docs/api/datasets.rst
@@ -0,0 +1,18 @@
+Datasets
+========
+
+Built-in dataset loaders for examples and testing. Each function returns a
+pandas ``Series`` with an appropriate index. See :doc:`/user-guide/datasets`
+for a summary table.
+
+.. autofunction:: dynaris.datasets.data.load_nile
+
+.. autofunction:: dynaris.datasets.data.load_airline
+
+.. autofunction:: dynaris.datasets.data.load_lynx
+
+.. autofunction:: dynaris.datasets.data.load_sunspots
+
+.. autofunction:: dynaris.datasets.data.load_temperature
+
+.. autofunction:: dynaris.datasets.data.load_gdp
diff --git a/docs/api/dlm.rst b/docs/api/dlm.rst
index 7da2f3d..64a7160 100644
--- a/docs/api/dlm.rst
+++ b/docs/api/dlm.rst
@@ -3,7 +3,16 @@ DLM --- High-Level Interface
 
 The :class:`~dynaris.dlm.api.DLM` class is the primary user-facing interface.
 It wraps a :class:`~dynaris.core.state_space.StateSpaceModel` with convenient
-``fit``, ``forecast``, ``smooth``, ``plot``, and ``summary`` methods.
+methods for the full modeling workflow:
+
+1. ``fit(y)`` --- run the Kalman filter
+2. ``smooth()`` --- run the RTS smoother
+3. ``forecast(steps)`` --- multi-step-ahead predictions
+4. ``plot(kind)`` --- visualize results
+5. ``summary()`` --- print model and fit information
+
+Most users only need this class. The lower-level filter, smoother, and
+forecast functions are available for advanced use cases.
 
 .. autoclass:: dynaris.dlm.api.DLM
    :members:
diff --git a/docs/api/estimation.rst b/docs/api/estimation.rst
index 4c9e86f..3cbfda1 100644
--- a/docs/api/estimation.rst
+++ b/docs/api/estimation.rst
@@ -1,11 +1,17 @@
 Parameter Estimation
 ====================
 
-Maximum likelihood estimation, EM algorithm, and model diagnostics.
+Maximum likelihood estimation, EM algorithm, residual diagnostics, and
+parameter transforms. See :doc:`/user-guide/estimation` for a guide on
+choosing between MLE and EM.
 
 MLE
 ---
 
+Gradient-based optimization of the log-likelihood using JAX autodiff.
+Flexible: supports any differentiable parameterization via a user-defined
+``model_fn``.
+
 .. autofunction:: dynaris.estimation.mle.fit_mle
 
 .. autoclass:: dynaris.estimation.mle.MLEResult
@@ -14,6 +20,9 @@ MLE
 EM Algorithm
 ------------
 
+Iterative variance estimation with guaranteed non-decreasing log-likelihood.
+Simpler setup than MLE --- just pass an initial model.
+
 .. autofunction:: dynaris.estimation.em.fit_em
 
 .. autoclass:: dynaris.estimation.em.EMResult
@@ -22,6 +31,8 @@ EM Algorithm
 Diagnostics
 -----------
 
+Tools for checking model adequacy after fitting.
+
 .. autofunction:: dynaris.estimation.diagnostics.standardized_residuals
 
 .. autofunction:: dynaris.estimation.diagnostics.acf
@@ -33,6 +44,8 @@ Diagnostics
 Transforms
 ----------
 
+Map unconstrained parameters to positive values for variance estimation.
+
 .. autofunction:: dynaris.estimation.transforms.softplus
 
 .. autofunction:: dynaris.estimation.transforms.inverse_softplus
diff --git a/docs/api/filters.rst b/docs/api/filters.rst
index 9f71e31..ed7bccf 100644
--- a/docs/api/filters.rst
+++ b/docs/api/filters.rst
@@ -1,7 +1,16 @@
 Kalman Filter
 =============
 
-Forward filtering for linear-Gaussian state-space models.
+Forward filtering for linear-Gaussian state-space models. The Kalman filter
+processes observations sequentially, computing the posterior state distribution
+at each time step.
+
+.. note::
+
+   Most users do not need to call these functions directly ---
+   :meth:`DLM.fit() <dynaris.dlm.api.DLM.fit>` wraps the Kalman filter
+   internally. These are available for advanced use cases requiring direct
+   access to intermediate filter quantities.
 
 .. autoclass:: dynaris.filters.kalman.KalmanFilter
    :members:
diff --git a/docs/api/forecast.rst b/docs/api/forecast.rst
index 41fb6b8..ab10f18 100644
--- a/docs/api/forecast.rst
+++ b/docs/api/forecast.rst
@@ -1,8 +1,17 @@
 Forecasting
 ===========
 
-Multi-step-ahead predictions with uncertainty quantification and
-batch processing.
+Multi-step-ahead predictions with uncertainty quantification and batch
+processing. Forecasts can be initialized from either the filtered or
+smoothed terminal state:
+
+- **From filtered state** (``forecast_from_filter``): uses observations up to
+  time :math:`T` only
+- **From smoothed state** (``forecast_from_smoother``): uses the full dataset
+  for a more refined starting point
+
+The high-level ``forecast`` function is called internally by
+:meth:`DLM.forecast() <dynaris.dlm.api.DLM.forecast>`.
 
 .. autofunction:: dynaris.forecast.forecast.forecast
 
@@ -18,3 +27,4 @@ batch processing.
 
 .. autoclass:: dynaris.forecast.forecast.ForecastResult
    :members:
+   :no-index:
diff --git a/docs/api/index.rst b/docs/api/index.rst
index 550f275..1e616e4 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -1,8 +1,33 @@
 API Reference
 =============
 
+Complete reference for all public classes and functions in dynaris.
+
++------------------+-------------------------------------------------------------+
+| Module           | Description                                                 |
++==================+=============================================================+
+| :doc:`dlm`       | High-level ``DLM`` class (fit, smooth, forecast, plot)      |
++------------------+-------------------------------------------------------------+
+| :doc:`components` | Six composable building blocks (``LocalLevel``, etc.)      |
++------------------+-------------------------------------------------------------+
+| :doc:`core`      | ``StateSpaceModel``, ``GaussianState``, result containers   |
++------------------+-------------------------------------------------------------+
+| :doc:`filters`   | Kalman filter (predict, update, full forward pass)          |
++------------------+-------------------------------------------------------------+
+| :doc:`smoothers` | Rauch-Tung-Striebel backward smoother                       |
++------------------+-------------------------------------------------------------+
+| :doc:`estimation` | MLE, EM algorithm, diagnostics, transforms                 |
++------------------+-------------------------------------------------------------+
+| :doc:`forecast`  | Multi-step forecasting and batch processing                 |
++------------------+-------------------------------------------------------------+
+| :doc:`plotting`  | Visualization functions for all plot kinds                  |
++------------------+-------------------------------------------------------------+
+| :doc:`datasets`  | Built-in dataset loaders                                    |
++------------------+-------------------------------------------------------------+
+
 .. toctree::
    :maxdepth: 2
+   :hidden:
 
    dlm
    components
@@ -12,3 +37,4 @@ API Reference
    estimation
    forecast
    plotting
+   datasets
diff --git a/docs/api/plotting.rst b/docs/api/plotting.rst
index 3cc3382..e9bd5c8 100644
--- a/docs/api/plotting.rst
+++ b/docs/api/plotting.rst
@@ -1,7 +1,15 @@
 Plotting
 ========
 
-Minimalist visualization functions using the cividis colormap.
+Visualization functions with a clean, minimalist style. All functions are
+called internally by :meth:`DLM.plot() <dynaris.dlm.api.DLM.plot>` but can
+also be used standalone.
+
+- ``plot_filtered`` --- observed data with filtered state overlay
+- ``plot_smoothed`` --- smoothed state estimates with confidence bands
+- ``plot_components`` --- decomposition into individual state components
+- ``plot_forecast`` --- forecast fan chart with historical context
+- ``plot_diagnostics`` --- QQ-plot, ACF, and histogram of residuals
 
 .. autofunction:: dynaris.plotting.plots.plot_filtered
 
diff --git a/docs/api/smoothers.rst b/docs/api/smoothers.rst
index 9e1006d..b37e47e 100644
--- a/docs/api/smoothers.rst
+++ b/docs/api/smoothers.rst
@@ -1,7 +1,15 @@
 RTS Smoother
 ============
 
-Rauch--Tung--Striebel backward smoother for linear-Gaussian models.
+Rauch-Tung-Striebel backward smoother for linear-Gaussian models. Uses the
+full dataset to refine state estimates, producing lower-variance posteriors
+than the forward-only Kalman filter.
+
+.. note::
+
+   Most users do not need to call these functions directly ---
+   :meth:`DLM.smooth() <dynaris.dlm.api.DLM.smooth>` wraps the RTS smoother
+   internally.
 
 .. autoclass:: dynaris.smoothers.rts.RTSSmoother
    :members:
diff --git a/docs/conf.py b/docs/conf.py
index 131c605..be1d57f 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -22,6 +22,7 @@
     "sphinx.ext.mathjax",
     "sphinx.ext.viewcode",
     "sphinx.ext.intersphinx",
+    "sphinx_copybutton",
 ]
 
 templates_path = ["_templates"]
@@ -29,18 +30,16 @@
 
 # -- Options for HTML output -------------------------------------------------
 
-html_theme = "alabaster"
+html_theme = "furo"
+html_title = "dynaris"
+html_logo = "_static/logo.png"
 html_theme_options = {
-    "description": "JAX-powered Dynamic Linear Models",
-    "github_user": "quant-sci",
-    "github_repo": "dynaris",
-    "github_button": True,
-    "fixed_sidebar": True,
-    "sidebar_collapse": True,
-    "page_width": "940px",
-    "sidebar_width": "220px",
+    "source_repository": "https://github.com/quant-sci/dynaris",
+    "source_branch": "main",
+    "source_directory": "docs/",
 }
 html_static_path = ["_static"]
+html_css_files = ["custom.css"]
 
 # -- Extension configuration -------------------------------------------------
 
diff --git a/docs/examples/index.rst b/docs/examples/index.rst
new file mode 100644
index 0000000..2adea4e
--- /dev/null
+++ b/docs/examples/index.rst
@@ -0,0 +1,55 @@
+Examples
+========
+
+Ready-to-run scripts demonstrating dynaris on real-world datasets. Each
+example lives in the ``examples/`` directory of the
+`repository <https://github.com/quant-sci/dynaris/tree/main/examples>`_.
+
+Classic / Introductory
+----------------------
+
+**Nile River Flow** --- ``nile_river.py``
+   Local level model on annual Nile discharge (1871--1970). The simplest
+   possible DLM: a random walk observed with noise. Demonstrates fitting,
+   smoothing, and forecasting.
+
+**Trend + Seasonality** --- ``trend_seasonal.py``
+   Simulated sales data with a local linear trend and monthly seasonality.
+   Shows how to compose components and extract individual contributions.
+
+Economics & Finance
+-------------------
+
+**Airline Passengers** --- ``airline_passengers.py``
+   Box-Jenkins airline data (1949--1960). Trend plus seasonal decomposition,
+   forecasting, and component visualization.
+
+**GDP Business Cycle** --- ``gdp_business_cycle.py``
+   US quarterly GDP growth. Local level plus AR(2) to capture business
+   cycle dynamics.
+
+Natural Sciences
+----------------
+
+**Lynx Population** --- ``lynx_population.py``
+   Canadian lynx trappings (1821--1934). Uses an autoregressive component to
+   model the ~10-year population cycle.
+
+**Sunspot Cycles** --- ``sunspot_cycles.py``
+   Annual sunspot numbers (1700--1987). Level plus cycle component to detect
+   the ~11-year solar cycle.
+
+**Global Temperature** --- ``global_temperature.py``
+   Annual temperature anomaly (1880--2023). Linear trend detection in climate
+   data.
+
+Advanced
+--------
+
+**Dynamic Regression** --- ``dynamic_regression.py``
+   Time-varying regression coefficients. Shows how regression weights evolve
+   over time when ``sigma_reg > 0``.
+
+**Panel Overview** --- ``panel_overview.py``
+   Demonstrates the multi-panel plot combining filtered, smoothed, forecast,
+   and diagnostic views in a single figure.
diff --git a/docs/getting-started/concepts.rst b/docs/getting-started/concepts.rst
new file mode 100644
index 0000000..6b14c12
--- /dev/null
+++ b/docs/getting-started/concepts.rst
@@ -0,0 +1,100 @@
+Key Concepts
+============
+
+What is a Dynamic Linear Model?
+--------------------------------
+
+A Dynamic Linear Model (DLM) is a state-space model where both the state
+evolution and observation equations are linear and driven by Gaussian noise.
+DLMs are a flexible framework for decomposing a time series into interpretable
+components --- trend, seasonality, regression effects, cycles --- and
+forecasting with uncertainty.
+
+At each time step :math:`t`, a DLM is defined by:
+
+- A **state** :math:`\boldsymbol{\theta}_t` that evolves over time (e.g., the
+  current level, slope, and seasonal effects)
+- An **observation** :math:`Y_t` that is a noisy linear function of the state
+
+The Kalman filter estimates the state given observations, and the RTS smoother
+refines those estimates using the full dataset.
+
+Components and composition
+--------------------------
+
+dynaris provides six building blocks, each defining a specific structure
+for the state-space matrices:
+
++-------------------------+------------+--------------------------------------------+
+| Component               | State dim  | Description                                |
++=========================+============+============================================+
+| ``LocalLevel``          | 1          | Random walk plus noise                     |
++-------------------------+------------+--------------------------------------------+
+| ``LocalLinearTrend``    | 2          | Level plus slope                           |
++-------------------------+------------+--------------------------------------------+
+| ``Seasonal``            | period - 1 | Dummy or Fourier seasonal effects          |
++-------------------------+------------+--------------------------------------------+
+| ``Regression``          | n_regressors | Dynamic or static coefficients           |
++-------------------------+------------+--------------------------------------------+
+| ``Autoregressive``      | order      | AR(p) in companion form                    |
++-------------------------+------------+--------------------------------------------+
+| ``Cycle``               | 2          | Damped stochastic sinusoid                 |
++-------------------------+------------+--------------------------------------------+
+
+Combine any of these with the ``+`` operator:
+
+.. code-block:: python
+
+   from dynaris import LocalLinearTrend, Seasonal, Cycle
+
+   model = LocalLinearTrend() + Seasonal(period=12) + Cycle(period=40, damping=0.95)
+   # Combined state_dim = 2 + 11 + 2 = 15
+
+Under the hood, this builds a single ``StateSpaceModel`` with block-diagonal
+system matrices (the **superposition principle** from West and Harrison, 1997).
+
+The workflow
+------------
+
+A typical dynaris session follows four steps:
+
+1. **Build** --- compose components into a model
+2. **Fit** --- run the Kalman filter on observed data
+3. **Smooth** --- run the RTS smoother for refined estimates
+4. **Forecast** --- project forward with uncertainty
+
+.. code-block:: python
+
+   from dynaris import LocalLevel, DLM
+
+   model = LocalLevel()
+   dlm = DLM(model)
+   dlm.fit(y)           # Kalman filter
+   dlm.smooth()         # RTS smoother
+   fc = dlm.forecast(steps=12)
+   dlm.plot(kind="panel")
+
+Notation
+--------
+
+dynaris follows the West and Harrison (1997) notation throughout:
+
++----------+-------------------------------+------------------------------------------+
+| Symbol   | Code                          | Meaning                                  |
++==========+===============================+==========================================+
+| **F**    | ``model.F`` / ``observation_matrix``  | Observation (regression) vector  |
++----------+-------------------------------+------------------------------------------+
+| **G**    | ``model.G`` / ``system_matrix``       | System (evolution) matrix        |
++----------+-------------------------------+------------------------------------------+
+| **V**    | ``model.V`` / ``obs_cov``             | Observational variance           |
++----------+-------------------------------+------------------------------------------+
+| **W**    | ``model.W`` / ``evolution_cov``       | Evolution covariance             |
++----------+-------------------------------+------------------------------------------+
+
+For the full mathematical treatment, see :doc:`/math`.
+
+References
+----------
+
+- West, M. and Harrison, J. (1997). *Bayesian Forecasting and Dynamic Models*,
+  2nd edition. Springer.
diff --git a/docs/getting-started/installation.rst b/docs/getting-started/installation.rst
new file mode 100644
index 0000000..b3179ca
--- /dev/null
+++ b/docs/getting-started/installation.rst
@@ -0,0 +1,48 @@
+Installation
+============
+
+From PyPI
+---------
+
+.. code-block:: bash
+
+   pip install dynaris
+
+Or with `uv <https://docs.astral.sh/uv/>`_:
+
+.. code-block:: bash
+
+   uv add dynaris
+
+From source
+-----------
+
+.. code-block:: bash
+
+   git clone https://github.com/quant-sci/dynaris.git
+   cd dynaris
+   pip install -e .
+
+Verify the installation
+-----------------------
+
+.. code-block:: python
+
+   import dynaris
+   print(dynaris.__version__)
+
+Dependencies
+------------
+
+dynaris requires Python 3.12+ and depends on:
+
+- `JAX <https://jax.readthedocs.io/>`_ (automatic differentiation and JIT)
+- NumPy
+- pandas
+- SciPy
+- Matplotlib
+
+.. note::
+
+   By default, JAX installs with CPU support. For GPU acceleration, see the
+   `JAX installation guide <https://jax.readthedocs.io/en/latest/installation.html>`_.
diff --git a/docs/getting-started/quickstart.rst b/docs/getting-started/quickstart.rst
new file mode 100644
index 0000000..5f298da
--- /dev/null
+++ b/docs/getting-started/quickstart.rst
@@ -0,0 +1,104 @@
+Quickstart
+==========
+
+This page walks through the core dynaris workflow: build a model, fit it to
+data, smooth, forecast, and plot.
+
+Your first DLM
+--------------
+
+The simplest DLM is a **local level** model --- a random walk observed with
+noise. Let's fit one to the classic Nile river dataset:
+
+.. code-block:: python
+
+   from dynaris import LocalLevel, DLM
+   from dynaris.datasets import load_nile
+
+   # 1. Load data (100 annual observations)
+   y = load_nile()
+
+   # 2. Define a local-level model
+   model = LocalLevel(sigma_level=38.0, sigma_obs=123.0)
+
+   # 3. Wrap in DLM and fit
+   dlm = DLM(model)
+   dlm.fit(y)
+
+   # 4. Inspect results
+   print(dlm.summary())
+
+The ``fit`` method runs the Kalman filter forward through the observations,
+computing filtered state estimates and the log-likelihood.
+
+Composing components
+--------------------
+
+The real power of dynaris is composition. Combine a trend with seasonality
+in a single line:
+
+.. code-block:: python
+
+   from dynaris import LocalLinearTrend, Seasonal, DLM
+   from dynaris.datasets import load_airline
+
+   y = load_airline()  # 144 monthly observations
+
+   model = LocalLinearTrend(sigma_level=2.0, sigma_slope=0.1) + Seasonal(period=12)
+   dlm = DLM(model)
+   dlm.fit(y)
+
+This produces a single state-space model with block-diagonal system matrices.
+See :doc:`/user-guide/components` for all six available components.
+
+Smoothing
+---------
+
+The Kalman filter processes observations forward in time. The RTS smoother
+uses future observations to refine past estimates:
+
+.. code-block:: python
+
+   dlm.fit(y).smooth()
+
+   # Smoothed states have lower variance than filtered states
+   df = dlm.smoothed_states_df()
+
+Forecasting
+-----------
+
+Generate multi-step-ahead forecasts with uncertainty intervals:
+
+.. code-block:: python
+
+   forecast_df = dlm.forecast(steps=24)
+   print(forecast_df)
+
+If you fit with a pandas Series that has a ``DatetimeIndex``, the forecast
+DataFrame continues the date index automatically.
+
+See :doc:`/user-guide/forecasting` for advanced options.
+
+Plotting
+--------
+
+View results with a single call:
+
+.. code-block:: python
+
+   # Single-figure overview (filtered, smoothed, forecast, diagnostics)
+   dlm.plot(kind="panel")
+
+   # Or individual plot types
+   dlm.plot(kind="filtered")
+   dlm.plot(kind="forecast", n_history=36)
+
+See :doc:`/user-guide/plotting` for all available plot kinds.
+
+Next steps
+----------
+
+- :doc:`concepts` --- understand the DLM framework and notation
+- :doc:`/user-guide/components` --- explore all six building blocks
+- :doc:`/user-guide/estimation` --- learn about parameter estimation (MLE and EM)
+- :doc:`/math` --- full mathematical foundations
diff --git a/docs/index.rst b/docs/index.rst
index 27db1cd..9c8cbbc 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,23 +1,101 @@
 dynaris
 =======
 
-A fast, composable, JAX-native DLM library with autodiff
-parameter estimation, multi-step forecasting, and clean visualization.
+**dynaris** is a JAX-powered Python library for Dynamic Linear Models (DLMs).
+Build composable state-space models, estimate parameters with automatic
+differentiation, and produce forecasts with uncertainty --- all in a few lines
+of code.
+
+.. code-block:: bash
+
+   pip install dynaris
 
 .. code-block:: python
 
-   from dynaris import LocalLevel, Seasonal, DLM
+   from dynaris import LocalLevel, DLM
+   from dynaris.datasets import load_nile
+
+   y = load_nile()
+   dlm = DLM(LocalLevel(sigma_level=38.0, sigma_obs=123.0))
+   dlm.fit(y).smooth()
+   dlm.forecast(steps=10)
+   dlm.plot(kind="panel")
+
+Why State-Space Models?
+-----------------------
+
+State-space models are a unifying framework for time series analysis. Rather
+than fitting a single equation to observed data, they maintain a latent
+**state** --- the true level, trend, seasonal pattern, or regression
+coefficient --- and update it recursively as new observations arrive. This
+separation of *what we observe* from *what drives the process* gives
+state-space models several structural advantages:
+
+- **Decomposition.** A complex series is expressed as the sum of interpretable
+  components (trend, seasonality, cycles, regression effects), each evolving
+  according to its own dynamics. The components are estimated jointly, not
+  sequentially stripped away.
+- **Exact uncertainty quantification.** The Kalman filter propagates a full
+  Gaussian posterior at every time step, so filtered estimates, smoothed
+  retrospectives, and multi-step forecasts all carry principled confidence
+  intervals --- no bootstrap or asymptotic approximation required.
+- **Missing data and irregular spacing.** When an observation is absent, the
+  filter simply skips the update step and lets the prior covariance grow.
+  No imputation heuristics are needed.
+- **Online and retrospective inference.** The same model supports real-time
+  filtering (estimate the present), smoothing (revise the past given the
+  future), and forecasting (project forward), all from a single set of
+  sufficient statistics.
 
-   model = LocalLevel() + Seasonal(period=12)
-   dlm = DLM(model)
-   dlm.fit(y)
-   dlm.forecast(steps=12)
-   dlm.plot()
+Dynamic Linear Models (DLMs) are the linear-Gaussian specialization of this
+framework, where closed-form Kalman recursions replace approximate inference.
+They are the workhorse behind structural time series decomposition, adaptive
+forecasting, and dynamic regression in fields from econometrics and signal
+processing to environmental monitoring.
+
+About dynaris
+-------------
+
+dynaris implements the full DLM inference pipeline in JAX. Models are built by
+composing six structural components --- ``LocalLevel``, ``LocalLinearTrend``,
+``Seasonal``, ``Cycle``, ``Regression``, and ``Autoregressive`` --- using the
+``+`` operator, which constructs block-diagonal state-space matrices via the
+superposition principle of West and Harrison (1997).
+
+Filtering and smoothing run inside ``jax.lax.scan``, making them JIT-compilable
+and end-to-end differentiable. This means the full Kalman log-likelihood can
+be optimized with gradient-based methods (``fit_mle``) or the EM algorithm
+(``fit_em``), and batch inference over many series parallelizes naturally
+through ``jax.vmap``.
 
 .. toctree::
    :maxdepth: 2
-   :caption: Contents
+   :caption: Getting Started
+
+   getting-started/installation
+   getting-started/quickstart
+   getting-started/concepts
+
+.. toctree::
+   :maxdepth: 2
+   :caption: User Guide
+
+   user-guide/index
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Theory
 
-   quickstart
    math
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Examples
+
+   examples/index
+
+.. toctree::
+   :maxdepth: 2
+   :caption: API Reference
+
    api/index
diff --git a/docs/math.rst b/docs/math.rst
index cca5854..de9eb65 100644
--- a/docs/math.rst
+++ b/docs/math.rst
@@ -5,6 +5,8 @@ This section presents the mathematical foundations of Dynamic Linear Models
 following the notation and framework of West and Harrison (1997),
 *Bayesian Forecasting and Dynamic Models*.
 
+For a plain-language introduction, see :doc:`getting-started/concepts`.
+
 The DLM Quadruple
 -----------------
 
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
deleted file mode 100644
index 19f0689..0000000
--- a/docs/quickstart.rst
+++ /dev/null
@@ -1,180 +0,0 @@
-Quickstart
-==========
-
-Installation
-------------
-
-.. code-block:: bash
-
-   pip install dynaris
-
-Your first DLM
---------------
-
-Dynaris models are built by composing components with the ``+`` operator,
-then wrapped in a :class:`~dynaris.dlm.api.DLM` for a high-level interface.
-
-.. code-block:: python
-
-   from dynaris import LocalLevel, DLM
-
-   # 1. Define a local-level model (random walk + noise)
-   model = LocalLevel(sigma_level=38.0, sigma_obs=123.0)
-
-   # 2. Wrap in DLM and fit
-   dlm = DLM(model)
-   dlm.fit(y)  # y can be a numpy array, JAX array, or pandas Series
-
-   # 3. Inspect results
-   print(dlm.summary())
-
-Composing components
---------------------
-
-The real power of dynaris is composition. Combine a trend with seasonality
-in a single line:
-
-.. code-block:: python
-
-   from dynaris import LocalLinearTrend, Seasonal, DLM
-
-   model = LocalLinearTrend(sigma_level=2.0, sigma_slope=0.1) + Seasonal(period=12)
-   dlm = DLM(model)
-   dlm.fit(y)
-
-This produces a single state-space model with block-diagonal system matrices.
-The state vector concatenates the trend states (level, slope) with the
-seasonal states (period - 1 dimensions).
-
-Smoothing
----------
-
-The Kalman filter processes observations forward in time. The RTS smoother
-uses future observations to refine past estimates:
-
-.. code-block:: python
-
-   dlm.fit(y).smooth()
-
-   # Smoothed states have lower variance than filtered states
-   df = dlm.smoothed_states_df()
-
-Forecasting
------------
-
-Generate multi-step-ahead forecasts with uncertainty intervals:
-
-.. code-block:: python
-
-   forecast_df = dlm.forecast(steps=24)
-   print(forecast_df)
-   #            mean     lower_95    upper_95
-   # 0     850.123     612.456    1087.790
-   # 1     850.123     598.234    1102.012
-   # ...
-
-If you fit with a pandas Series that has a ``DatetimeIndex``, the forecast
-DataFrame continues the date index automatically.
-
-Plotting
---------
-
-All plots use a clean, minimalist style with the cividis colormap:
-
-.. code-block:: python
-
-   # Filtered vs observed
-   dlm.plot(kind="filtered")
-
-   # Smoothed states
-   dlm.plot(kind="smoothed")
-
-   # Forecast fan chart
-   dlm.forecast(steps=24)
-   dlm.plot(kind="forecast", n_history=36)
-
-   # Residual diagnostics (QQ-plot, ACF, histogram)
-   dlm.plot(kind="diagnostics")
-
-   # Component decomposition (requires smooth first)
-   dlm.smooth()
-   dlm.plot(kind="components", component_dims={
-       "Level": 0,
-       "Slope": 1,
-       "Seasonal": 2,
-   })
-
-Parameter estimation
---------------------
-
-Estimate unknown variance parameters via maximum likelihood:
-
-.. code-block:: python
-
-   import jax.numpy as jnp
-   from dynaris import LocalLevel
-   from dynaris.estimation import fit_mle
-
-   def model_fn(params):
-       return LocalLevel(
-           sigma_level=jnp.exp(params[0]),
-           sigma_obs=jnp.exp(params[1]),
-       )
-
-   result = fit_mle(model_fn, y, init_params=jnp.zeros(2))
-   print(f"Optimal log-likelihood: {result.log_likelihood:.2f}")
-   fitted_model = result.model
-
-Or use the EM algorithm:
-
-.. code-block:: python
-
-   from dynaris.estimation import fit_em
-
-   result = fit_em(y, initial_model, max_iter=100)
-   print(f"Converged: {result.converged}")
-
-Batch processing
-----------------
-
-Fit or forecast multiple series in parallel with ``jax.vmap``:
-
-.. code-block:: python
-
-   import jax.numpy as jnp
-
-   # y_batch: shape (n_series, T, obs_dim)
-   batch_result = dlm.fit_batch(y_batch)
-   print(batch_result.log_likelihood)  # shape (n_series,)
-
-All components
---------------
-
-Dynaris provides six composable DLM building blocks:
-
-+-------------------------+------------+--------------------------------------------+
-| Component               | State dim  | Description                                |
-+=========================+============+============================================+
-| ``LocalLevel``          | 1          | Random walk + noise                        |
-+-------------------------+------------+--------------------------------------------+
-| ``LocalLinearTrend``    | 2          | Level + slope                              |
-+-------------------------+------------+--------------------------------------------+
-| ``Seasonal``            | period - 1 | Dummy or Fourier form                      |
-+-------------------------+------------+--------------------------------------------+
-| ``Regression``          | n_regressors | Dynamic/static coefficients              |
-+-------------------------+------------+--------------------------------------------+
-| ``Autoregressive``      | order      | AR(p) in companion form                    |
-+-------------------------+------------+--------------------------------------------+
-| ``Cycle``               | 2          | Damped stochastic sinusoid                 |
-+-------------------------+------------+--------------------------------------------+
-
-Combine any of these with ``+``:
-
-.. code-block:: python
-
-   model = (
-       LocalLinearTrend()
-       + Seasonal(period=12)
-       + Cycle(period=40, damping=0.95)
-   )
-   # state_dim = 2 + 11 + 2 = 15
diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 0000000..f442e2a
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,3 @@
+sphinx>=9.1.0
+furo>=2024.8.6
+sphinx-copybutton>=0.5
diff --git a/docs/user-guide/batch-processing.rst b/docs/user-guide/batch-processing.rst
new file mode 100644
index 0000000..dc2a93f
--- /dev/null
+++ b/docs/user-guide/batch-processing.rst
@@ -0,0 +1,50 @@
+Batch Processing
+================
+
+dynaris supports fitting and forecasting multiple time series in parallel
+using ``jax.vmap``.
+
+Batch fitting
+-------------
+
+Pass a batch of series as a 3D array with shape ``(n_series, T, obs_dim)``:
+
+.. code-block:: python
+
+   import jax.numpy as jnp
+   from dynaris import LocalLevel, DLM
+
+   model = LocalLevel()
+   dlm = DLM(model)
+
+   # y_batch: (n_series, T, 1)
+   batch_result = dlm.fit_batch(y_batch)
+   print(batch_result.log_likelihood)  # shape (n_series,)
+
+Each series is filtered independently, but all series run in parallel on the
+same hardware (CPU cores or GPU).
+
+Batch forecasting
+-----------------
+
+After batch fitting, generate forecasts for all series at once:
+
+.. code-block:: python
+
+   from dynaris.forecast import forecast_batch
+
+   fc = forecast_batch(batch_result, model, steps=12)
+
+Low-level API
+-------------
+
+The batch functions wrap ``jax.vmap`` over the single-series equivalents:
+
+.. code-block:: python
+
+   from dynaris.forecast import fit_batch, forecast_batch
+
+   batch_filter = fit_batch(model, y_batch)
+   batch_fc = forecast_batch(batch_filter, model, steps=12)
+
+See :doc:`/api/forecast` for the full API.
diff --git a/docs/user-guide/components.rst b/docs/user-guide/components.rst
new file mode 100644
index 0000000..af56e54
--- /dev/null
+++ b/docs/user-guide/components.rst
@@ -0,0 +1,161 @@
+Components
+==========
+
+dynaris provides six composable DLM building blocks. Each function returns a
+:class:`~dynaris.core.state_space.StateSpaceModel` that can be combined with
+others using the ``+`` operator.
+
++-------------------------+------------+--------------------------------------------+
+| Component               | State dim  | Description                                |
++=========================+============+============================================+
+| ``LocalLevel``          | 1          | Random walk plus noise                     |
++-------------------------+------------+--------------------------------------------+
+| ``LocalLinearTrend``    | 2          | Level plus slope                           |
++-------------------------+------------+--------------------------------------------+
+| ``Seasonal``            | period - 1 | Dummy or Fourier seasonal effects          |
++-------------------------+------------+--------------------------------------------+
+| ``Regression``          | n_regressors | Dynamic or static coefficients           |
++-------------------------+------------+--------------------------------------------+
+| ``Autoregressive``      | order      | AR(p) in companion form                    |
++-------------------------+------------+--------------------------------------------+
+| ``Cycle``               | 2          | Damped stochastic sinusoid                 |
++-------------------------+------------+--------------------------------------------+
+
+Trend components
+----------------
+
+LocalLevel
+~~~~~~~~~~
+
+The simplest DLM: a random walk observed with noise.
+
+.. math::
+
+   \mu_t = \mu_{t-1} + \omega_t, \quad
+   Y_t = \mu_t + \nu_t
+
+.. code-block:: python
+
+   from dynaris import LocalLevel
+
+   model = LocalLevel(sigma_level=1.0, sigma_obs=1.0)
+
+**When to use:** stationary or slowly changing series without a clear trend
+or seasonality. Classic example: Nile river annual flow.
+
+LocalLinearTrend
+~~~~~~~~~~~~~~~~
+
+Extends ``LocalLevel`` with a slope (growth rate) component.
+
+.. math::
+
+   \mu_t = \mu_{t-1} + \beta_{t-1} + \omega_{\mu,t}, \quad
+   \beta_t = \beta_{t-1} + \omega_{\beta,t}
+
+.. code-block:: python
+
+   from dynaris import LocalLinearTrend
+
+   model = LocalLinearTrend(sigma_level=2.0, sigma_slope=0.1, sigma_obs=1.0)
+
+**When to use:** series with a changing trend direction, such as GDP growth
+or temperature anomalies.
+
+Periodic components
+-------------------
+
+Seasonal
+~~~~~~~~
+
+Models repeating patterns at a fixed period. Supports both **dummy** (default)
+and **Fourier** forms.
+
+.. code-block:: python
+
+   from dynaris import Seasonal
+
+   # Dummy-form seasonal (default)
+   model = Seasonal(period=12, sigma_seasonal=0.5)
+
+   # Fourier-form seasonal
+   model = Seasonal(period=12, sigma_seasonal=0.5, form="fourier")
+
+**When to use:** monthly, quarterly, or weekly data with repeating patterns.
+State dimension is ``period - 1``.
+
+See :doc:`/math` for the mathematical formulation.
+
+Cycle
+~~~~~
+
+A damped stochastic sinusoid for quasi-periodic behavior.
+
+.. code-block:: python
+
+   from dynaris import Cycle
+
+   model = Cycle(period=40, damping=0.95, sigma_cycle=1.0)
+
+**When to use:** series with approximate cycles of known period, such as
+sunspot activity (~11 years) or business cycles. Setting ``damping=1.0``
+gives an undamped cycle; values below 1.0 let the cycle decay.
+
+Other components
+----------------
+
+Regression
+~~~~~~~~~~
+
+Dynamic (time-varying) or static regression coefficients.
+
+.. code-block:: python
+
+   from dynaris import Regression
+
+   # Dynamic coefficients (random walk)
+   model = Regression(n_regressors=2, sigma_reg=0.1)
+
+   # Static coefficients (set sigma_reg=0)
+   model = Regression(n_regressors=2, sigma_reg=0.0)
+
+**When to use:** when external predictors (regressors) influence the series.
+With ``sigma_reg > 0``, coefficients evolve over time; with ``sigma_reg = 0``,
+they are constant.
+
+Autoregressive
+~~~~~~~~~~~~~~
+
+An AR(p) process in companion-form state space.
+
+.. code-block:: python
+
+   from dynaris import Autoregressive
+
+   model = Autoregressive(coefficients=[0.5, -0.3], sigma_ar=1.0)
+
+**When to use:** series with serial correlation not captured by trend or
+seasonal components. Common in population dynamics (e.g., lynx data).
+
+Composing components
+--------------------
+
+Combine any components with ``+`` to build richer models:
+
+.. code-block:: python
+
+   from dynaris import LocalLinearTrend, Seasonal, Cycle, DLM
+
+   model = (
+       LocalLinearTrend(sigma_level=2.0, sigma_slope=0.1)
+       + Seasonal(period=12, sigma_seasonal=0.5)
+       + Cycle(period=40, damping=0.95)
+   )
+   # state_dim = 2 + 11 + 2 = 15
+
+   dlm = DLM(model)
+   dlm.fit(y)
+
+This uses the DLM **superposition principle**: system matrices become
+block-diagonal, and observation noise adds across components. See
+:doc:`/math` for details.
diff --git a/docs/user-guide/datasets.rst b/docs/user-guide/datasets.rst
new file mode 100644
index 0000000..84952de
--- /dev/null
+++ b/docs/user-guide/datasets.rst
@@ -0,0 +1,40 @@
+Built-in Datasets
+=================
+
+dynaris ships with six classic time series datasets for examples and testing.
+Each loader returns a pandas ``Series`` with an appropriate index.
+
++------------------+--------------------+------+-----------+-----------------+
+| Dataset          | Loader             | N    | Frequency | Domain          |
++==================+====================+======+===========+=================+
+| Nile river flow  | ``load_nile()``    | 100  | Annual    | Hydrology       |
++------------------+--------------------+------+-----------+-----------------+
+| Airline          | ``load_airline()`` | 144  | Monthly   | Transportation  |
+| passengers       |                    |      |           |                 |
++------------------+--------------------+------+-----------+-----------------+
+| Lynx population  | ``load_lynx()``    | 114  | Annual    | Ecology         |
++------------------+--------------------+------+-----------+-----------------+
+| Sunspot numbers  | ``load_sunspots()``| 288  | Annual    | Astronomy       |
++------------------+--------------------+------+-----------+-----------------+
+| Global           | ``load_temperature | 144  | Annual    | Climate         |
+| temperature      | ()``               |      |           |                 |
++------------------+--------------------+------+-----------+-----------------+
+| US GDP growth    | ``load_gdp()``     | 319  | Quarterly | Economics       |
++------------------+--------------------+------+-----------+-----------------+
+
+Usage
+-----
+
+.. code-block:: python
+
+   from dynaris.datasets import load_airline
+
+   y = load_airline()
+   print(y.head())
+   print(f"Shape: {y.shape}, Index: {y.index[0]} to {y.index[-1]}")
+
+All loaders accept no arguments and return a pandas ``Series``. The index
+is a ``DatetimeIndex`` for monthly/quarterly data or an integer index for
+annual data.
+
+See :doc:`/api/datasets` for the full API reference.
diff --git a/docs/user-guide/estimation.rst b/docs/user-guide/estimation.rst
new file mode 100644
index 0000000..487bbd9
--- /dev/null
+++ b/docs/user-guide/estimation.rst
@@ -0,0 +1,118 @@
+Parameter Estimation
+====================
+
+dynaris provides two approaches for estimating unknown variance parameters:
+**maximum likelihood estimation (MLE)** via automatic differentiation and the
+**EM algorithm**.
+
+Maximum Likelihood (MLE)
+------------------------
+
+The log-likelihood is computed by the Kalman filter's prediction error
+decomposition. Since the entire computation runs in JAX, gradients are
+obtained via ``jax.grad`` and passed to a gradient-based optimizer.
+
+.. code-block:: python
+
+   import jax.numpy as jnp
+   from dynaris import LocalLevel
+   from dynaris.estimation import fit_mle
+
+   def model_fn(params):
+       return LocalLevel(
+           sigma_level=jnp.exp(params[0]),
+           sigma_obs=jnp.exp(params[1]),
+       )
+
+   result = fit_mle(model_fn, y, init_params=jnp.zeros(2))
+   print(f"Log-likelihood: {result.log_likelihood:.2f}")
+   fitted_model = result.model
+
+The ``model_fn`` maps unconstrained parameters to a ``StateSpaceModel``.
+Use ``jnp.exp`` (log transform) or ``softplus`` to ensure variance parameters
+stay positive.
+
+The result is an :class:`~dynaris.estimation.mle.MLEResult` containing the
+optimized model, final parameters, and log-likelihood.
+
+EM Algorithm
+------------
+
+The Expectation-Maximization algorithm alternates between running the Kalman
+smoother (E-step) and updating variance estimates (M-step):
+
+.. code-block:: python
+
+   from dynaris.estimation import fit_em
+
+   result = fit_em(y, initial_model, max_iter=100)
+   print(f"Converged: {result.converged}")
+   print(f"Iterations: {result.n_iter}")
+   fitted_model = result.model
+
+The result is an :class:`~dynaris.estimation.em.EMResult` with the fitted
+model, convergence status, and iteration count.
+
+MLE vs EM
+---------
+
++-------------------+----------------------------------+----------------------------------+
+| Criterion         | MLE                              | EM                               |
++===================+==================================+==================================+
+| Speed             | Fewer iterations (gradient info) | More iterations, but each is     |
+|                   |                                  | cheap and closed-form            |
++-------------------+----------------------------------+----------------------------------+
+| Flexibility       | Any differentiable               | Variance parameters only         |
+|                   | parameterization                 |                                  |
++-------------------+----------------------------------+----------------------------------+
+| Convergence       | Can find local minima            | Guaranteed non-decreasing        |
+|                   |                                  | log-likelihood                   |
++-------------------+----------------------------------+----------------------------------+
+| Setup             | Requires writing ``model_fn``    | Just pass the initial model      |
++-------------------+----------------------------------+----------------------------------+
+
+**Rule of thumb:** use MLE for complex parameterizations; use EM when you only
+need variance estimates and want a simple interface.
+
+Parameter transforms
+--------------------
+
+Variance parameters must be positive. dynaris provides two transforms for
+mapping unconstrained parameters to valid ranges:
+
+- **Log transform:** :math:`\sigma^2 = \exp(\psi)` --- simple, widely used
+- **Softplus transform:** :math:`\sigma^2 = \log(1 + \exp(\psi))` --- smoother
+  gradient near zero
+
+.. code-block:: python
+
+   from dynaris.estimation import softplus, inverse_softplus
+
+   # Map unconstrained -> positive
+   sigma_sq = softplus(raw_param)
+
+   # Map positive -> unconstrained
+   raw_param = inverse_softplus(sigma_sq)
+
+Diagnostics
+-----------
+
+After fitting, check model adequacy with residual diagnostics:
+
+.. code-block:: python
+
+   from dynaris.estimation import standardized_residuals, acf, ljung_box
+
+   resid = standardized_residuals(filter_result, model)
+   autocorr = acf(resid, max_lag=20)
+   lb = ljung_box(resid, max_lag=10)
+   print(f"Ljung-Box p-value: {lb.p_value:.4f}")
+
+Or use the built-in diagnostic plot:
+
+.. code-block:: python
+
+   dlm.plot(kind="diagnostics")
+
+This produces a QQ-plot, ACF plot, and histogram of standardized residuals.
+See :doc:`/api/estimation` for the full API.
diff --git a/docs/user-guide/forecasting.rst b/docs/user-guide/forecasting.rst
new file mode 100644
index 0000000..55983c5
--- /dev/null
+++ b/docs/user-guide/forecasting.rst
@@ -0,0 +1,78 @@
+Forecasting
+===========
+
+dynaris produces multi-step-ahead forecasts with uncertainty intervals by
+iterating the state-space prior equations forward without new observations.
+
+Basic forecasting
+-----------------
+
+After fitting (and optionally smoothing), call ``forecast``:
+
+.. code-block:: python
+
+   from dynaris import LocalLinearTrend, Seasonal, DLM
+   from dynaris.datasets import load_airline
+
+   y = load_airline()
+   model = LocalLinearTrend() + Seasonal(period=12)
+   dlm = DLM(model)
+   dlm.fit(y).smooth()
+
+   forecast_df = dlm.forecast(steps=24)
+   print(forecast_df)
+   #            mean     lower_95    upper_95
+   # ...
+
+The returned DataFrame contains the forecast mean and 95% confidence bands.
+
+DatetimeIndex propagation
+-------------------------
+
+If you fit with a pandas Series that has a ``DatetimeIndex``, the forecast
+DataFrame continues the date index:
+
+.. code-block:: python
+
+   # y has monthly DatetimeIndex 1949-01 to 1960-12
+   forecast_df = dlm.forecast(steps=12)
+   # Index continues: 1961-01, 1961-02, ...
+
+Filtered vs smoothed initialization
+------------------------------------
+
+Forecasts can start from either the filtered or smoothed terminal state:
+
+- **Filtered** (default after ``fit``): uses only past observations up to
+  time :math:`T`
+- **Smoothed** (after ``smooth``): uses the full dataset, giving a more
+  refined starting point
+
+The lower-level functions give explicit control:
+
+.. code-block:: python
+
+   from dynaris.forecast import forecast_from_filter, forecast_from_smoother
+
+   fc_filt = forecast_from_filter(filter_result, model, steps=12)
+   fc_smooth = forecast_from_smoother(smoother_result, model, steps=12)
+
+Confidence bands
+----------------
+
+Confidence bands widen with the forecast horizon as evolution noise
+accumulates:
+
+.. code-block:: python
+
+   from dynaris.forecast import confidence_bands
+
+   lower, upper = confidence_bands(forecast_result, level=0.95)
+
+Visualize with:
+
+.. code-block:: python
+
+   dlm.plot(kind="forecast", n_history=36)
+
+See :doc:`/api/forecast` for the full API.
diff --git a/docs/user-guide/index.rst b/docs/user-guide/index.rst
new file mode 100644
index 0000000..4a74a32
--- /dev/null
+++ b/docs/user-guide/index.rst
@@ -0,0 +1,15 @@
+User Guide
+==========
+
+This guide covers each part of dynaris in detail --- from model components
+and parameter estimation to forecasting, plotting, and batch processing.
+
+.. toctree::
+   :maxdepth: 2
+
+   components
+   estimation
+   forecasting
+   plotting
+   batch-processing
+   datasets
diff --git a/docs/user-guide/plotting.rst b/docs/user-guide/plotting.rst
new file mode 100644
index 0000000..1f18186
--- /dev/null
+++ b/docs/user-guide/plotting.rst
@@ -0,0 +1,99 @@
+Plotting
+========
+
+dynaris provides minimalist visualization functions. All plots are accessible
+through the ``DLM.plot()`` method via the ``kind`` parameter.
+
+Plot kinds
+----------
+
+filtered
+~~~~~~~~
+
+Overlay filtered state estimates on the observed data.
+
+.. code-block:: python
+
+   dlm.fit(y)
+   dlm.plot(kind="filtered")
+
+Shows the Kalman filter's one-step-ahead estimates with confidence intervals.
+
+smoothed
+~~~~~~~~
+
+Display smoothed (retrospective) state estimates.
+
+.. code-block:: python
+
+   dlm.fit(y).smooth()
+   dlm.plot(kind="smoothed")
+
+Smoothed estimates have lower variance because they use the full dataset.
+
+forecast
+~~~~~~~~
+
+Fan chart showing the forecast mean and confidence bands, with recent
+historical observations for context.
+
+.. code-block:: python
+
+   dlm.forecast(steps=24)
+   dlm.plot(kind="forecast", n_history=36)
+
+The ``n_history`` parameter controls how many past observations appear.
+
+diagnostics
+~~~~~~~~~~~
+
+Residual diagnostic panel with QQ-plot, ACF, and histogram.
+
+.. code-block:: python
+
+   dlm.plot(kind="diagnostics")
+
+Use this to check whether the model's assumptions hold.
+
+components
+~~~~~~~~~~
+
+Decompose the series into individual state components. Requires smoothing
+first and a mapping of component names to state dimensions:
+
+.. code-block:: python
+
+   dlm.smooth()
+   dlm.plot(kind="components", component_dims={
+       "Level": 0,
+       "Slope": 1,
+       "Seasonal": 2,
+   })
+
+panel
+~~~~~
+
+A single-figure overview combining filtered, smoothed, forecast, and
+diagnostics:
+
+.. code-block:: python
+
+   dlm.fit(y).smooth()
+   dlm.forecast(steps=12)
+   dlm.plot(kind="panel")
+
+Customization
+-------------
+
+All plot methods accept an optional ``ax`` parameter to draw on an existing
+Matplotlib axes, and a ``title`` parameter:
+
+.. code-block:: python
+
+   import matplotlib.pyplot as plt
+
+   fig, ax = plt.subplots()
+   dlm.plot(kind="filtered", ax=ax, title="Nile River Flow")
+   plt.show()
+
+See :doc:`/api/plotting` for the full function signatures.
diff --git a/pyproject.toml b/pyproject.toml
index 3baa10b..17da597 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,8 @@ dev = [
     "mypy>=1.14",
     "pandas-stubs>=3.0.0.260204",
     "sphinx>=9.1.0",
-    "alabaster>=1.0.0",
+    "furo>=2024.8.6",
+    "sphinx-copybutton>=0.5",
 ]
 
 [build-system]
diff --git a/uv.lock b/uv.lock
index 0085df1..902e2ee 100644
--- a/uv.lock
+++ b/uv.lock
@@ -13,6 +13,18 @@ resolution-markers = [
     "python_full_version < '3.13' and sys_platform != 'emscripten' and sys_platform != 'win32'",
 ]
 
+[[package]]
+name = "accessible-pygments"
+version = "0.0.5"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "pygments" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/bc/c1/bbac6a50d02774f91572938964c582fff4270eee73ab822a4aeea4d8b11b/accessible_pygments-0.0.5.tar.gz", hash = "sha256:40918d3e6a2b619ad424cb91e556bd3bd8865443d9f22f1dcdf79e33c8046872", size = 1377899, upload-time = "2024-05-10T11:23:10.216Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/8d/3f/95338030883d8c8b91223b4e21744b04d11b161a3ef117295d8241f50ab4/accessible_pygments-0.0.5-py3-none-any.whl", hash = "sha256:88ae3211e68a1d0b011504b2ffc1691feafce124b845bd072ab6f9f66f34d4b7", size = 1395903, upload-time = "2024-05-10T11:23:08.421Z" },
+]
+
 [[package]]
 name = "alabaster"
 version = "1.0.0"
@@ -31,6 +43,19 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/77/f5/21d2de20e8b8b0408f0681956ca2c69f1320a3848ac50e6e7f39c6159675/babel-2.18.0-py3-none-any.whl", hash = "sha256:e2b422b277c2b9a9630c1d7903c2a00d0830c409c59ac8cae9081c92f1aeba35", size = 10196845, upload-time = "2026-02-01T12:30:53.445Z" },
 ]
 
+[[package]]
+name = "beautifulsoup4"
+version = "4.14.3"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "soupsieve" },
+    { name = "typing-extensions" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/c3/b0/1c6a16426d389813b48d95e26898aff79abbde42ad353958ad95cc8c9b21/beautifulsoup4-4.14.3.tar.gz", hash = "sha256:6292b1c5186d356bba669ef9f7f051757099565ad9ada5dd630bd9de5fa7fb86", size = 627737, upload-time = "2025-11-30T15:08:26.084Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/1a/39/47f9197bdd44df24d67ac8893641e16f386c984a0619ef2ee4c51fbbc019/beautifulsoup4-4.14.3-py3-none-any.whl", hash = "sha256:0918bfe44902e6ad8d57732ba310582e98da931428d231a5ecb9e7c703a735bb", size = 107721, upload-time = "2025-11-30T15:08:24.087Z" },
+]
+
 [[package]]
 name = "certifi"
 version = "2026.2.25"
@@ -292,7 +317,7 @@ wheels = [
 
 [[package]]
 name = "dynaris"
-version = "0.1.0"
+version = "0.1.1"
 source = { editable = "." }
 dependencies = [
     { name = "jax" },
@@ -305,13 +330,14 @@ dependencies = [
 
 [package.dev-dependencies]
 dev = [
-    { name = "alabaster" },
+    { name = "furo" },
     { name = "mypy" },
     { name = "pandas-stubs" },
     { name = "pytest" },
     { name = "pytest-cov" },
     { name = "ruff" },
     { name = "sphinx" },
+    { name = "sphinx-copybutton" },
 ]
 
 [package.metadata]
@@ -326,13 +352,14 @@ requires-dist = [
 
 [package.metadata.requires-dev]
 dev = [
-    { name = "alabaster", specifier = ">=1.0.0" },
+    { name = "furo", specifier = ">=2024.8.6" },
     { name = "mypy", specifier = ">=1.14" },
     { name = "pandas-stubs", specifier = ">=3.0.0.260204" },
     { name = "pytest", specifier = ">=8.0" },
     { name = "pytest-cov", specifier = ">=6.0" },
     { name = "ruff", specifier = ">=0.9" },
     { name = "sphinx", specifier = ">=9.1.0" },
+    { name = "sphinx-copybutton", specifier = ">=0.5" },
 ]
 
 [[package]]
@@ -376,6 +403,22 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/fd/ba/56147c165442cc5ba7e82ecf301c9a68353cede498185869e6e02b4c264f/fonttools-4.62.1-py3-none-any.whl", hash = "sha256:7487782e2113861f4ddcc07c3436450659e3caa5e470b27dc2177cade2d8e7fd", size = 1152647, upload-time = "2026-03-13T13:54:22.735Z" },
 ]
 
+[[package]]
+name = "furo"
+version = "2025.12.19"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "accessible-pygments" },
+    { name = "beautifulsoup4" },
+    { name = "pygments" },
+    { name = "sphinx" },
+    { name = "sphinx-basic-ng" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/ec/20/5f5ad4da6a5a27c80f2ed2ee9aee3f9e36c66e56e21c00fde467b2f8f88f/furo-2025.12.19.tar.gz", hash = "sha256:188d1f942037d8b37cd3985b955839fea62baa1730087dc29d157677c857e2a7", size = 1661473, upload-time = "2025-12-19T17:34:40.889Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/f4/b2/50e9b292b5cac13e9e81272c7171301abc753a60460d21505b606e15cf21/furo-2025.12.19-py3-none-any.whl", hash = "sha256:bb0ead5309f9500130665a26bee87693c41ce4dbdff864dbfb6b0dae4673d24f", size = 339262, upload-time = "2025-12-19T17:34:38.905Z" },
+]
+
 [[package]]
 name = "idna"
 version = "3.11"
@@ -1220,6 +1263,15 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/c8/78/3565d011c61f5a43488987ee32b6f3f656e7f107ac2782dd57bdd7d91d9a/snowballstemmer-3.0.1-py3-none-any.whl", hash = "sha256:6cd7b3897da8d6c9ffb968a6781fa6532dce9c3618a4b127d920dab764a19064", size = 103274, upload-time = "2025-05-09T16:34:50.371Z" },
 ]
 
+[[package]]
+name = "soupsieve"
+version = "2.8.3"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/7b/ae/2d9c981590ed9999a0d91755b47fc74f74de286b0f5cee14c9269041e6c4/soupsieve-2.8.3.tar.gz", hash = "sha256:3267f1eeea4251fb42728b6dfb746edc9acaffc4a45b27e19450b676586e8349", size = 118627, upload-time = "2026-01-20T04:27:02.457Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/46/2c/1462b1d0a634697ae9e55b3cecdcb64788e8b7d63f54d923fcd0bb140aed/soupsieve-2.8.3-py3-none-any.whl", hash = "sha256:ed64f2ba4eebeab06cc4962affce381647455978ffc1e36bb79a545b91f45a95", size = 37016, upload-time = "2026-01-20T04:27:01.012Z" },
+]
+
 [[package]]
 name = "sphinx"
 version = "9.1.0"
@@ -1248,6 +1300,30 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/73/f7/b1884cb3188ab181fc81fa00c266699dab600f927a964df02ec3d5d1916a/sphinx-9.1.0-py3-none-any.whl", hash = "sha256:c84fdd4e782504495fe4f2c0b3413d6c2bf388589bb352d439b2a3bb99991978", size = 3921742, upload-time = "2025-12-31T15:09:25.561Z" },
 ]
 
+[[package]]
+name = "sphinx-basic-ng"
+version = "1.0.0b2"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "sphinx" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/98/0b/a866924ded68efec7a1759587a4e478aec7559d8165fac8b2ad1c0e774d6/sphinx_basic_ng-1.0.0b2.tar.gz", hash = "sha256:9ec55a47c90c8c002b5960c57492ec3021f5193cb26cebc2dc4ea226848651c9", size = 20736, upload-time = "2023-07-08T18:40:54.166Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/3c/dd/018ce05c532a22007ac58d4f45232514cd9d6dd0ee1dc374e309db830983/sphinx_basic_ng-1.0.0b2-py3-none-any.whl", hash = "sha256:eb09aedbabfb650607e9b4b68c9d240b90b1e1be221d6ad71d61c52e29f7932b", size = 22496, upload-time = "2023-07-08T18:40:52.659Z" },
+]
+
+[[package]]
+name = "sphinx-copybutton"
+version = "0.5.2"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "sphinx" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/fc/2b/a964715e7f5295f77509e59309959f4125122d648f86b4fe7d70ca1d882c/sphinx-copybutton-0.5.2.tar.gz", hash = "sha256:4cf17c82fb9646d1bc9ca92ac280813a3b605d8c421225fd9913154103ee1fbd", size = 23039, upload-time = "2023-04-14T08:10:22.998Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/9e/48/1ea60e74949eecb12cdd6ac43987f9fd331156388dcc2319b45e2ebb81bf/sphinx_copybutton-0.5.2-py3-none-any.whl", hash = "sha256:fb543fd386d917746c9a2c50360c7905b605726b9355cd26e9974857afeae06e", size = 13343, upload-time = "2023-04-14T08:10:20.844Z" },
+]
+
 [[package]]
 name = "sphinxcontrib-applehelp"
 version = "2.0.0"