diff --git a/rdtools/test/analysis_chains_test.py b/rdtools/test/analysis_chains_test.py
index efe78df4..e4e48a68 100644
--- a/rdtools/test/analysis_chains_test.py
+++ b/rdtools/test/analysis_chains_test.py
@@ -87,6 +87,59 @@ def sensor_analysis_aggregated_no_filter(sensor_parameters):
     return rd_analysis
 
 
+@pytest.fixture
+def clearsky_example_data(basic_parameters):
+    # Import the example data
+    file_url = (
+        "https://datahub.duramat.org/dataset/"
+        "a49bb656-7b36-437a-8089-1870a40c2a7d/"
+        "resource/d2c3fcf4-4f5f-47ad-8743-fc29"
+        "f1356835/download/pvdaq_system_4_2010"
+        "-2016_subset_soil_signal.csv"
+    )
+    cache_file = "PVDAQ_system_4_2010-2016_subset_soilsignal.pickle"
+
+    try:
+        df = pd.read_pickle(cache_file)
+    except FileNotFoundError:
+        df = pd.read_csv(file_url, index_col=0, parse_dates=True)
+        df.to_pickle(cache_file)
+
+    # Specify the Metadata
+    meta = {
+        "latitude": 39.7406,
+        "longitude": -105.1774,
+        "timezone": "Etc/GMT+7",
+        "gamma_pdc": -0.005,
+        "azimuth": 180,
+        "tilt": 40,
+        "power_dc_rated": 1000.0,
+        "temp_model_params": "open_rack_glass_polymer",
+    }
+
+    # Set the timezone
+    df.index = df.index.tz_localize(meta["timezone"])
+
+    # Select two years of data
+    df_crop = df[df.index < (df.index[0] + pd.Timedelta(days=2 * 365 + 1))]
+
+    basic_parameters["pv"] = df_crop["ac_power"]
+    basic_parameters["poa_global"] = df_crop["poa_irradiance"]
+    basic_parameters["temperature_ambient"] = df_crop["ambient_temp"]
+    basic_parameters["interp_freq"] = "1min"
+
+    # Set the pvlib location
+    loc = pvlib.location.Location(meta["latitude"], meta["longitude"], tz=meta["timezone"])
+
+    cs_input = dict(
+        pvlib_location=loc,
+        pv_tilt=meta["tilt"],
+        pv_azimuth=meta["azimuth"],
+        solar_position_method="ephemeris",  # just to improve test execution speed
+    )
+    return basic_parameters, cs_input
+
+
 def test_interpolation(basic_parameters, degradation_trend):
 
     power = degradation_trend
@@ -436,9 +489,7 @@ def test_no_gamma_pdc(sensor_parameters):
 
 
 @pytest.fixture
-def clearsky_parameters(
-    basic_parameters, sensor_parameters, cs_input, degradation_trend
-):
+def clearsky_parameters(basic_parameters, sensor_parameters, cs_input, degradation_trend):
     # clear-sky weather data.  Uses TrendAnalysis's internal clear-sky
     # functions to generate the data.
     rd_analysis = TrendAnalysis(**sensor_parameters)
@@ -461,6 +512,16 @@ def clearsky_analysis(cs_input, clearsky_parameters):
     return rd_analysis
 
 
+@pytest.fixture
+def clearsky_pvlib_analysis(clearsky_example_data):
+    clearsky_parameters_example, cs_input_example = clearsky_example_data
+    rd_analysis = TrendAnalysis(**clearsky_parameters_example)
+    rd_analysis.set_clearsky(**cs_input_example)
+    rd_analysis.filter_params["clearsky_filter"] = {"model": "pvlib"}
+    rd_analysis.clearsky_analysis(analyses=["yoy_degradation"])
+    return rd_analysis
+
+
 @pytest.fixture
 def clearsky_optional(cs_input, clearsky_analysis):
     # optional parameters to exercise other branches
@@ -486,15 +547,33 @@ def sensor_clearsky_analysis(cs_input, clearsky_parameters):
     return rd_analysis
 
 
+@pytest.fixture
+def sensor_clearsky_pvlib_analysis(clearsky_example_data):
+    clearsky_parameters_example, cs_input_example = clearsky_example_data
+    rd_analysis = TrendAnalysis(**clearsky_parameters_example)
+    rd_analysis.set_clearsky(**cs_input_example)
+    rd_analysis.filter_params = {}  # disable all index-based filters
+    rd_analysis.filter_params["sensor_clearsky_filter"] = {"model": "pvlib"}
+    rd_analysis.sensor_analysis(analyses=["yoy_degradation"])
+    return rd_analysis
+
+
 def test_clearsky_analysis(clearsky_analysis):
     yoy_results = clearsky_analysis.results["clearsky"]["yoy_degradation"]
     ci = yoy_results["rd_confidence_interval"]
     rd = yoy_results["p50_rd"]
-    print(ci)
     assert pytest.approx(rd, abs=1e-2) == -5.15
     assert pytest.approx(ci, abs=1e-2) == [-5.17, -5.13]
 
 
+def test_clearsky_pvlib_analysis(clearsky_pvlib_analysis):
+    yoy_results = clearsky_pvlib_analysis.results["clearsky"]["yoy_degradation"]
+    ci = yoy_results["rd_confidence_interval"]
+    rd = yoy_results["p50_rd"]
+    assert pytest.approx(rd, abs=1e-2) == -1.589
+    assert pytest.approx(ci, abs=1e-2) == [-2.417, -0.861]
+
+
 def test_clearsky_analysis_filter_components(clearsky_analysis):
     columns = clearsky_analysis.clearsky_filter_components_aggregated.columns
     assert {'two_way_window_filter'} == set(columns)
@@ -523,11 +602,18 @@ def test_sensor_clearsky_analysis(sensor_clearsky_analysis):
     yoy_results = sensor_clearsky_analysis.results["sensor"]["yoy_degradation"]
     ci = yoy_results["rd_confidence_interval"]
     rd = yoy_results["p50_rd"]
-    print(ci)
     assert -5.18 == pytest.approx(rd, abs=1e-2)
     assert [-5.18, -5.18] == pytest.approx(ci, abs=1e-2)
 
 
+def test_sensor_clearsky_pvlib_analysis(sensor_clearsky_pvlib_analysis):
+    yoy_results = sensor_clearsky_pvlib_analysis.results["sensor"]["yoy_degradation"]
+    ci = yoy_results["rd_confidence_interval"]
+    rd = yoy_results["p50_rd"]
+    assert -1.478 == pytest.approx(rd, abs=1e-2)
+    assert [-1.926, -0.649] == pytest.approx(ci, abs=1e-2)
+
+
 @pytest.fixture
 def clearsky_analysis_exp_power(clearsky_parameters, clearsky_optional):
     power_expected = normalization.pvwatts_dc_power(