diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 39bc164f..e79e3b5b 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -36,7 +36,7 @@ jobs:
       - name: Build wheels
         run: |
           expected=0
-          for pyver in 39 310 311 312 313; do
+          for pyver in 39 310 311 312 313 314; do
             pybin="/opt/python/cp${pyver}-cp${pyver}/bin/python"
             if [ ! -f "$pybin" ]; then
               echo "ERROR: Expected Python interpreter not found: $pybin"
@@ -65,7 +65,7 @@ jobs:
     runs-on: macos-14
     strategy:
       matrix:
-        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13', '3.14']
     steps:
       - uses: actions/checkout@v4
 
@@ -95,7 +95,7 @@ jobs:
     runs-on: windows-latest
     strategy:
       matrix:
-        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13', '3.14']
     steps:
       - uses: actions/checkout@v4
 
diff --git a/.github/workflows/rust-test.yml b/.github/workflows/rust-test.yml
index d22c1d85..60c76e26 100644
--- a/.github/workflows/rust-test.yml
+++ b/.github/workflows/rust-test.yml
@@ -81,7 +81,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest, ubuntu-24.04-arm]
-        python-version: ['3.11', '3.13']
+        python-version: ['3.11', '3.13', '3.14']
 
     steps:
       - uses: actions/checkout@v4
@@ -184,7 +184,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: '3.13'
+          python-version: '3.14'
 
       - name: Install dependencies
         # Keep in sync with pyproject.toml [project.dependencies] and [project.optional-dependencies.dev]
diff --git a/README.md b/README.md
index 42cc866c..63793bd3 100644
--- a/README.md
+++ b/README.md
@@ -2769,7 +2769,7 @@ Returns DataFrame with columns: `unit`, `quality_score`, `outcome_trend_score`,
 
 ## Requirements
 
-- Python 3.9 - 3.13
+- Python 3.9 - 3.14
 - numpy >= 1.20
 - pandas >= 1.3
 - scipy >= 1.7
diff --git a/pyproject.toml b/pyproject.toml
index 192d0a96..38ac4704 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -8,7 +8,7 @@ version = "3.0.0"
 description = "Difference-in-Differences causal inference with sklearn-like API. Callaway-Sant'Anna, Synthetic DiD, Honest DiD, event studies, parallel trends."
 readme = "README.md"
 license = "MIT"
-requires-python = ">=3.9,<3.14"
+requires-python = ">=3.9,<3.15"
 authors = [
     {name = "diff-diff contributors"}
 ]
@@ -39,6 +39,7 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     "Topic :: Scientific/Engineering :: Mathematics",
     "Topic :: Scientific/Engineering :: Information Analysis",
     "Topic :: Scientific/Engineering",
diff --git a/rust/Cargo.toml b/rust/Cargo.toml
index 26ee6ace..67426d09 100644
--- a/rust/Cargo.toml
+++ b/rust/Cargo.toml
@@ -2,6 +2,7 @@
 name = "diff_diff_rust"
 version = "3.0.0"
 edition = "2021"
+rust-version = "1.84"
 description = "Rust backend for diff-diff DiD library"
 license = "MIT"
 
@@ -20,10 +21,10 @@ accelerate = ["ndarray/blas", "dep:blas-src", "blas-src/accelerate"]
 openblas = ["ndarray/blas"]
 
 [dependencies]
-# PyO3 0.22 supports Python 3.8-3.13
-pyo3 = "0.22"
-numpy = "0.22"
-ndarray = { version = "0.16", features = ["rayon"] }
+# PyO3 0.28 supports Python 3.9-3.14
+pyo3 = "0.28"
+numpy = "0.28"
+ndarray = { version = "0.17", features = ["rayon"] }
 rand = "0.8"
 rand_xoshiro = "0.6"
 rayon = "1.8"
diff --git a/rust/src/bootstrap.rs b/rust/src/bootstrap.rs
index e6539286..71d5093f 100644
--- a/rust/src/bootstrap.rs
+++ b/rust/src/bootstrap.rs
@@ -48,7 +48,7 @@ pub fn generate_bootstrap_weights_batch<'py>(
         }
     };
 
-    Ok(weights.to_pyarray_bound(py))
+    Ok(weights.to_pyarray(py))
 }
 
 /// Generate Rademacher weights: ±1 with equal probability.
diff --git a/rust/src/linalg.rs b/rust/src/linalg.rs
index 1260fa6b..8d4c5854 100644
--- a/rust/src/linalg.rs
+++ b/rust/src/linalg.rs
@@ -144,20 +144,20 @@ pub fn solve_ols<'py>(
             // Rank-deficient: cannot compute valid vcov, return NaN matrix
             let mut nan_vcov = Array2::<f64>::zeros((k, k));
             nan_vcov.fill(f64::NAN);
-            Some(nan_vcov.to_pyarray_bound(py))
+            Some(nan_vcov.to_pyarray(py))
         } else {
             // Full rank: compute robust vcov normally
             let cluster_arr = cluster_ids.as_ref().map(|c| c.as_array().to_owned());
             let vcov_arr = compute_robust_vcov_internal(&x_arr, &residuals.view(), cluster_arr.as_ref(), n, k)?;
-            Some(vcov_arr.to_pyarray_bound(py))
+            Some(vcov_arr.to_pyarray(py))
         }
     } else {
         None
     };
 
     Ok((
-        coefficients.to_pyarray_bound(py),
-        residuals.to_pyarray_bound(py),
+        coefficients.to_pyarray(py),
+        residuals.to_pyarray(py),
         vcov,
     ))
 }
@@ -186,7 +186,7 @@ pub fn compute_robust_vcov<'py>(
     let n = x_arr.nrows();
     let k = x_arr.ncols();
     let vcov = compute_robust_vcov_internal(&x_arr, &residuals_arr, cluster_arr.as_ref(), n, k)?;
-    Ok(vcov.to_pyarray_bound(py))
+    Ok(vcov.to_pyarray(py))
 }
 
 /// Internal implementation of robust variance-covariance computation.
diff --git a/rust/src/trop.rs b/rust/src/trop.rs
index 605303a7..4257d4a8 100644
--- a/rust/src/trop.rs
+++ b/rust/src/trop.rs
@@ -44,7 +44,7 @@ pub fn compute_unit_distance_matrix<'py>(
 
     let dist_matrix = compute_unit_distance_matrix_internal(&y_arr, &d_arr);
 
-    Ok(dist_matrix.to_pyarray_bound(py))
+    Ok(dist_matrix.to_pyarray(py))
 }
 
 /// Internal implementation of unit distance matrix computation.
@@ -1098,7 +1098,7 @@ pub fn bootstrap_trop_variance<'py>(
     };
 
     let estimates_arr = Array1::from_vec(bootstrap_estimates);
-    Ok((estimates_arr.to_pyarray_bound(py), se))
+    Ok((estimates_arr.to_pyarray(py), se))
 }
 
 // ============================================================================
@@ -1838,7 +1838,7 @@ pub fn bootstrap_trop_variance_global<'py>(
     };
 
     let estimates_arr = Array1::from_vec(bootstrap_estimates);
-    Ok((estimates_arr.to_pyarray_bound(py), se))
+    Ok((estimates_arr.to_pyarray(py), se))
 }
 
 #[cfg(test)]
diff --git a/rust/src/weights.rs b/rust/src/weights.rs
index 001b9b23..00c95b01 100644
--- a/rust/src/weights.rs
+++ b/rust/src/weights.rs
@@ -54,7 +54,7 @@ pub fn compute_synthetic_weights<'py>(
     let weights =
         compute_synthetic_weights_internal(&y_control_arr, &y_treated_arr, lambda_reg, max_iter, tol)?;
 
-    Ok(weights.to_pyarray_bound(py))
+    Ok(weights.to_pyarray(py))
 }
 
 /// Internal implementation of synthetic weight computation.
@@ -137,7 +137,7 @@ pub fn project_simplex<'py>(
 ) -> PyResult<Bound<'py, PyArray1<f64>>> {
     let v_arr = v.as_array();
     let result = project_simplex_internal(&v_arr);
-    Ok(result.to_pyarray_bound(py))
+    Ok(result.to_pyarray(py))
 }
 
 /// Internal implementation of simplex projection.
@@ -607,7 +607,7 @@ pub fn sc_weight_fw<'py>(
         min_decrease,
         max_iter,
     );
-    Ok(result.to_pyarray_bound(py))
+    Ok(result.to_pyarray(py))
 }
 
 /// Compute SDID time weights via Frank-Wolfe optimization.
@@ -637,7 +637,7 @@ pub fn compute_time_weights<'py>(
     let y_post = y_post_control.as_array();
 
     let result = compute_time_weights_internal(&y_pre, &y_post, zeta_lambda, intercept, min_decrease, max_iter_pre_sparsify, max_iter);
-    Ok(result.to_pyarray_bound(py))
+    Ok(result.to_pyarray(py))
 }
 
 pub(crate) fn compute_time_weights_internal(
@@ -720,7 +720,7 @@ pub fn compute_sdid_unit_weights<'py>(
         &y_pre, &y_tr_mean, zeta_omega, intercept, min_decrease,
         max_iter_pre_sparsify, max_iter,
     );
-    Ok(result.to_pyarray_bound(py))
+    Ok(result.to_pyarray(py))
 }
 
 pub(crate) fn compute_sdid_unit_weights_internal(