Movatterモバイル変換

Type:	Package
Title:	Time-Varying Minimum Variance Portfolio
Version:	1.0.5
Date:	2025-06-27
Description:	Provides the estimation of a time-dependent covariance matrix of returns with the intended use for portfolio optimization. The package offers methods for determining the optimal number of factors to be used in the covariance estimation, a hypothesis test of time-varying covariance, and user-friendly functions for portfolio optimization and rolling window evaluation. The local PCA method, method for determining the number of factors, and associated hypothesis test are based on Su and Wang (2017) <doi:10.1016/j.jeconom.2016.12.004>. The approach to time-varying portfolio optimization follows Fan et al. (2024) <doi:10.1016/j.jeconom.2022.08.007>. The regularisation applied to the residual covariance matrix adopts the technique introduced by Chen et al. (2019) <doi:10.1016/j.jeconom.2019.04.025>.
License:	MIT + file LICENSE
URL:	https://github.com/erilill/TV-MVP
BugReports:	https://github.com/erilill/TV-MVP/issues
Encoding:	UTF-8
Depends:	R (≥ 3.6.0)
Imports:	R6, cli, prettyunits, dplyr, ggplot2, tidyr
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0)
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2025-06-27 15:24:41 UTC; erikl_xzy542i
Author:	Erik Lillrank [aut, cre], Yukai Yang [aut]
Maintainer:	Erik Lillrank <erik.lillrank@gmail.com>
Repository:	CRAN
Date/Publication:	2025-06-27 15:50:06 UTC

TVMVP: Time-Varying Minimum Variance Portfolio Optimization

Description

The TVMVP package provides tools for estimating time-dependent covariance matrices using kernel-weighted principal component analysis. These estimates can then be used for portfolio optimization in a dynamic setting.

Details

The method involves five steps: (1) determining the number of factors, (2) estimating kernel-weighted PCA, (3) regularizing the idiosyncratic error covariance, (4) estimating the total covariance matrix, and (5) computing optimal portfolio weights.

An optional step includes a hypothesis test to check whether the covariance matrix is time-invariant.

The local PCA method, method for determining the number of factors, and associated hypothesis test are based on Su and Wang (2017). The approach to time-varying portfolio optimization follows Fan et al. (2024). The regularisation applied to the residual covariance matrix adopts the technique introduced by Chen et al. (2019).

The methodology implemented in this package closely follows Fan et al. (2024).The original authors provide a Matlab implementation at https://github.com/RuikeWu/TV-MVP.

The default kernel function in the package is the Epanechnikov kernel. Otherkernel functions can also be used, however these are not implemented in thepackage. In order to do this, write an R function with an integrablekernel function and use this as input in the functions with argumentkernel_func. It should be constructed ascustom_kernel <- function(u){...}.

Similarly, the bandwidth function which is implemented in the package is the Silverman's rule of thumb. For most functions, simply setbandwidthto your preferred bandwidth, however forexpanding_tvmvp,only Silverman's is implemented in this version of the package.

Authors and Maintainer

Authors: Erik Lillrank and Yukai Yang
Maintainer: Erik Lillrank
Department of Statistics, Uppsala University
erik.lillrank@gmail.com,yukai.yang@statistik.uu.se

Functions

determine_factors

Selects the optimal number of factors via an information criterion.

hyptest

Hypothesis test for time-invariant covariance matrices. Bootstrap p-values supported.

predict_portfolio

Optimizes portfolio weights for out-of-sample prediction of portfolio performance.

expanding_tvmvp

Evaluates MVP performance in a expanding window framework.

time_varying_cov

Estimates the time-varying covariance matrix.

silverman

Silverman's rule of thumb bandwidth formula.

Class

TVMVP

Time Varying Minimum Variance Portfolio (TVMVP) Class.

References

Lillrank, E. (2025).A Time-Varying Factor Approach to Covariance Estimation

Su, L., & Wang, X. (2017). On time-varying factor models: Estimation and testing. Journal of Econometrics, 198(1), 84–101.

Fan, Q., Wu, R., Yang, Y., & Zhong, W. (2024). Time-varying minimum variance portfolio. Journal of Econometrics, 239(2), 105339.

Chen, J., Li, D., & Linton, O. (2019). A new semiparametric estimation approach for large dynamic covariance matrices with multiple conditioning variables. Journal of Econometrics, 212(1), 155–176.

Author(s)

Maintainer: Erik Lillrankerik.lillrank@gmail.com (ORCID)

Authors:

• Yukai Yangyukai.yang@statistik.uu.se (ORCID)

See Also

Useful links:

• https://github.com/erilill/TV-MVP

• Report bugs athttps://github.com/erilill/TV-MVP/issues

Time Varying Minimum Variance Portfolio (TVMVP) Class

Description

This class implements a time-varying minimum variance portfolio using locallysmoothed principal component analysis (PCA) to estimate the time-dependentcovariance matrix.

This class provides a flexible interface to:

• Set return data ($set_data())

• Determine optimal number of factors ($determine_factors())

• Conduct test of constant factor loadings ($hyptest())

• Time-dependent covariance estimation ($time_varying_cov())

• Portfolio optimization ($predict_portfolio())

• Expanding window evaluation ($expanding_tvmvp())

• Extract cached results ($get_optimal_m(),$get_IC_values(),$get_bootstrap())

Looking for package description? SeeTVMVP-package.

Usage

# Initial object of class TVMVPtv <- TVMVP$new()# Set datatv$set_data(returns) # Returns must be T times p matrix# Determine number of factorstv$determine_factors(max_m=10)# Test for constant loadingstv$hyptest()# Estimate time-dependent covariance matrixcov <- tv$time_varying_cov()# Evaluate TVMVP performance on historical datamvp_results <- tv$expanding_tvmvp(                 initial_window = 60,                 rebal_period   = 5,                 max_factors    = 10,                 return_type    = "daily")                 # Make out-of-sample prediction and compute weightspredictions <- tv$predict_portfolio(horizon=5,                                   min_return = 0.01,                                   max_SR = TRUE)# Extract weightspredictions$getWeights("MVP")

#'

Arguments

data

T × p (time periods by assets) matrix of returns.

bandwidth

Numerical. Bandwidth parameter used for local smoothing in the local PCA

max_m

Integer. Maximum number of factors to be tested when determining the optimal number of factors.

optimal_m

The optimal number of factors to use in covariance estimation.

Methods

$new(data = NULL)

Initialize object of class TVMVP. Optionally pass returns matrix.

$set_data(data)

Set the data. Must beT × p (time periods by assets) matrix.

$get_data()

Get the data.

$set()

Manually set arguments of the object.

$determine_factors()

Determines optimal number of factors based on BIC-type information criterion.

$get_optimal_m{}

Prints optimal number of factors,optimal_m.

$get_IC_values()

Prints IC-values for the different number of factors tested usingdetermine_factors.

$hyptest()

Hypothesis test of constant loadings.

$get_bootstrap()

Prints bootstrap test statistics from the hypothesis test.

$predict_portfolio()

Optimizes portfolio weights for out-of-sample prediction of portfolio performance.

$expanding_tvmvp()

Evaluates MVP performance in a expanding window framework.

$time_varying_cov()

Estimates the time-varying covariance matrix.

$silverman()

Silverman's rule of thumb bandwidth formula.

Adaptive Selection of the Shrinkage Parameter\rho for POET

Description

This function selects an optimal shrinkage parameter\rho for the residual covarianceestimation procedure. It does so by dividing the data into groups and comparing a shrunk covariancematrix (computed on one subsample) to a benchmark covariance (computed on another subsample) usingthe Frobenius norm. The candidate\rho that minimizes the total squared Frobenius norm differenceis selected.

Usage

adaptive_poet_rho(  R,  M0 = 10,  rho_grid = seq(0.001, 2, length.out = 20),  epsilon2 = 1e-06)

Arguments

Details

The function proceeds as follows:

1. The total number of observationsT is halved to defineT_1 andT_2. Specifically:

   T_1 = \left\lfloor \frac{T}{2} \times \left(1 - \frac{1}{\log(T)}\right) \right\rfloor

   T_2 = \left\lfloor \frac{T}{2} \right\rfloor - T_1

2. The sample is divided into\left\lfloor T/(2M_0) \right\rfloor groups (withM_0 observations left out in between).

3. For each group, two subsamples are defined:

   • Subsample 1: the firstT_1 observations of the group.

   • Subsample 2: the lastT_2 observations of the group, after skippingM_0 observations following subsample 1.

4. For each group and a given candidate\rho inrho_grid, the covariance matrixS_1is computed from subsample 1, and then shrunk using soft-thresholding:

   S_{1,\text{shrunk}} = \text{soft\_threshold}\left(S_1, \rho \times \text{mean}\left(|S_1|_{\text{off-diag}}\right)\right)

5. The total squared Frobenius norm betweenS_{1,\text{shrunk}}and the covariance matrixS_2 (from subsample 2) is computed across all groups.

6. The function scansrho_grid to find the\rho minimizing total error. Additionally, it computes\rho_1as\epsilon_2 plus the smallest\rho for which the smallest eigenvalue of the shrunk covariance is positive.

Value

A list containing:

• best_rho: The selected optimal shrinkage parameter\hat{\rho} that minimizes the totalsquared Frobenius norm difference.

• rho_1: The lower bound for\rho derived from the minimum eigenvalue criteria (adjusted byepsilon2).

• min_Fnorm: The minimum total squared Frobenius norm difference achieved.

Boundary Kernel Function

Description

This function computes boundary kernel weights for a given time periodt within a datasetof sizeT. It adjusts the kernel weights near the boundaries to account for edge effects,ensuring that the weights sum to one.

Usage

boundary_kernel(t, r, iT, h, kernel_func)

Arguments

Details

The boundary kernel function adjusts kernel weights near the start and end of the dataset to mitigateedge effects commonly encountered in kernel-based methods. The function performs the following steps:

1. Scales the difference between the current timet and reference timer by theproduct of total time periodsT and bandwidthh.

2. Applies the kernel function to the scaled difference and adjusts by the bandwidth.

3. Determines if the current time period is within the lower or upper boundary regions based onT_h = \lfloor T \times h \rfloor.

4. Computes the integral of the kernel function over the adjusted limits to ensure the weightssum to one in boundary regions.

Value

A numeric scalar representing the boundary-adjusted kernel weight for the given time period.

Function to compute expected returns using a simple model selection approach

Description

Function to compute expected returns using a simple model selection approach

Usage

comp_expected_returns(returns, horizon)

Arguments

ComputeB_{pT} Statistic for Covariance Time-Variation Hypothesis Testing

Description

This function calculates theB_{pT} statistic, which is part of the hypothesistesting procedure to determine whether the covariance matrix of asset returns is time-varying.It incorporates kernel-weighted local and global factor interactions along with residuals.

Usage

compute_B_pT(local_factors, global_factors, residuals, h, iT, ip, kernel_func)

Arguments

Details

The function performs the following steps:

1. Computes the sum of squared residuals for each time periods.

2. Constructs the kernel matrixK[s,t] by applying theboundary_kernelfunction to each pair of time periods(s,t).

3. Calculates the local dot-product matrixL[s,t] as the dot product betweenthe local factors at times andt.

4. Computes the global dot-product matrixG[s,t] as the dot product betweenthe global factors at times andt.

5. Computes the element-wise squared difference betweenK * L andG,multiplies it by the residuals, and sums over alls,t.

6. Scales the aggregated value by\frac{\sqrt{h}}{T^2 \sqrt{p}} to obtainB_{pT}.

Value

A numeric scalarB_{pT} representing the computed statistic based onkernel-weighted factor interactions and residuals.

ComputeM_{\hat{}} Statistic for Covariance Time-Variation Hypothesis Testing

Description

This function calculates theM_{\hat{}} statistic, which measures the average squareddiscrepancy between local and global factor models across all assets and time periods.It quantifies the difference between locally estimated factors/loadings and their globalcounterparts.

Usage

compute_M_hat(  local_factors,  global_factors,  local_loadings,  global_loadings,  iT,  ip,  m)

Arguments

Details

The function performs the following steps:

1. Initializes theM_{\hat{}} statistic to zero.

2. If the number of factorsm is equal to one, it ensures thatglobal_loadings andglobal_factors are treated as matrices.

3. Iterates over each asseti = 1 toN and each time periodt = 1 toT.

4. For each asset and time period, computes:

   • common_H1: The dot product of the local loadings and local factors.

   • common_H0: The dot product of the global loadings and global factors.

   • The squared difference(common\_H1 - common\_H0)^2 and adds it toM_{\hat{}}.

5. After all iterations, normalizesM_{\hat{}} by dividing by the product ofN andT.

Value

A numeric scalarM_{\hat{}} representing the average squared discrepancybetween local and global factor models across all assets and time periods.

ComputeV_{pT} Statistic for Covariance Time-Variation Hypothesis Testing

Description

This function calculates theV_{pT} statistic, which is part of the hypothesistesting procedure to determine whether the covariance matrix of asset returns is time-varying.It incorporates kernel-weighted factor interactions and residual correlations.

Usage

compute_V_pT(local_factors, residuals, h, iT, ip, kernel_func)

Arguments

Details

The function performs the following steps:

1. Iterates over each pair of time periods(s, r) wheres < r.

2. Computes the two-fold convolution kernel value\bar{K}_{sr} using thetwo_fold_convolution_kernel function.

3. Calculates the squared dot product of local factors weighted by the factor covariancematrix.

4. Computes the squared dot product of residuals between time periodss andr.

5. Aggregates these values across all relevant time period pairs and scales by\frac{2}{T^2 × p × h} to obtainV_{pT}.

Value

A numeric scalarV_{pT} representing the computed statistic based onkernel-weighted factor interactions and residual correlations.

Compute Sigma_0 p.93 Su and Wang (2017).

Description

Compute Sigma_0 p.93 Su and Wang (2017).

Usage

compute_sigma_0(res_set, iT, ip)

Arguments

Value

Sigma_0 from page 93 in Su and Wang (2017).

Determine the Optimal Number of Factors via an Information Criterion

Description

This function selects the optimal number of factors for a local principal componentanalysis (PCA) model of asset returns. It computes an BIC-type information criterion (IC) for each candidatenumber of factors, based on the sum of squared residuals (SSR) from the PCA reconstruction and apenalty term that increases with the number of factors. The optimal number of factors is chosen as theone that minimizes the IC. The procedure is available either as a stand-alonefunction or as a method in the 'TVMVP' R6 class.

Usage

determine_factors(returns, max_m, bandwidth = silverman(returns))

Arguments

Details

Two usage styles:

# Function interfacedetermine_factors(returns, max_m = 5)# R6 method interfacetv <- TVMVP$new()tv$set_data(returns)tv$determine_factors(max_m = 5)tv$get_optimal_m()tv$get_IC_values()

When using the method form, if 'max_m' or 'bandwidth' are omitted,they default to values stored in the object. Results are cached andretrievable via class methods.

For each candidate number of factorsm (from 1 tomax_m), the function:

1. Performs a local PCA on the returns at each time pointr = 1,\dots,T usingm factors.

2. Computes a reconstruction of the returns and the corresponding residuals:

   \text{Residual}_r = R_r - F_r \Lambda_r,

   whereR_r is the return at timer, andF_r and\Lambda_r are the local factors and loadings, respectively.

3. Computes the average sum of squared residuals (SSR) as:

   V(m) = \frac{1}{pT} \sum_{r=1}^{T} \| \text{Residual}_r \|^2.

4. Adds a penalty term that increases withR:

   \text{Penalty}(m) = m × \frac{(p + T × \text{bandwidth})}{(pT × \text{bandwidth})} \log\left(\frac{pT × \text{bandwidth}}{(p + T × \text{bandwidth})}\right).

5. The information criterion is defined as:

   \text{IC}(m) = \log\big(V(m)\big) + \text{Penalty}(m).

The optimal number of factors is then chosen as the value ofm that minimizes\text{IC}(m).

Value

A list with:

• optimal_m: Integer. The optimal number of factors.

• IC_values: Numeric vector of IC values for each candidatem.

References

Su, L., & Wang, X. (2017). On time-varying factor models: Estimation and testing. Journal of Econometrics, 198(1), 84–101.

Examples

set.seed(123)returns <- matrix(rnorm(100 * 30), nrow = 100, ncol = 30)# Function usageresult <- determine_factors(returns, max_m = 5)print(result$optimal_m)print(result$IC_values)# R6 usagetv <- TVMVP$new()tv$set_data(returns)tv$determine_factors(max_m = 5)tv$get_optimal_m()tv$get_IC_values()

Epanechnikov Kernel Function

Description

This function computes the value of the Epanechnikov kernel for a given inputu.The Epanechnikov kernel is a popular choice in kernel density estimation due to its optimalproperties in minimizing mean integrated squared error.

Usage

epanechnikov_kernel(u)

Arguments

Details

The Epanechnikov kernel is defined as:

K(u) = \begin{cases}\frac{3}{4}(1 - u^2) & \text{if } |u| \leq 1, \\0 & \text{otherwise}.\end{cases}

This function applies the above formula to each element of the input vectoru.

Value

A numeric vector of kernel values corresponding to each inputu.

Examples

# Define a range of u valuesu_values <- seq(-1.5, 1.5, by = 0.1)# Compute Epanechnikov kernel valueskernel_values <- epanechnikov_kernel(u_values)

Estimate Local Covariance

Description

This internal function computes a time-varying covariance matrix estimate for a givenwindow of asset returns by combining factor-based and sparse residual covariance estimation.It uses results from a local PCA to form residuals and then applies an adaptive thresholdingprocedure (viaadaptive_poet_rho()) to shrink the residual covariance.

Usage

estimate_residual_cov_poet_local(  localPCA_results,  returns,  M0 = 10,  rho_grid = seq(0.005, 2, length.out = 30),  floor_value = 1e-12,  epsilon2 = 1e-06)

Arguments

Details

The function follows these steps:

1. **Local Residuals:**Extract the local loadings\Lambda_t from the last element oflocalPCA_results\$loadings andfactors\hat{F} fromlocalPCA_results\$f_hat. Letw_t denote the corresponding kernel weights.The local residuals are computed as:

   U_{\text{local}} = R - F \Lambda_t,

   whereR is the returns matrix.

2. **Adaptive Thresholding:**The function callsadaptive_poet_rho() onU_{\text{local}}to select an optimal shrinkage parameter\hat{\rho}_t.

3. **Residual Covariance Estimation:**The raw residual covariance is computed as:

   S_{u,\text{raw}} = \frac{1}{T} U_{\text{local}}^\top U_{\text{local}},

   and a threshold is set as:

   \text{threshold} = \hat{\rho}_t × \text{mean}(|S_{u,\text{raw}}|),

   where the mean is taken over the off-diagonal elements.Soft-thresholding is then applied to obtain the shrunk residual covariance matrix\hat{S}_u.

4. **Total Covariance Estimation:**The final covariance matrix is constructed by combining the factor component with the shrunk residual covariance:

   \Sigma_R(t) = \Lambda_t \left(\frac{F^\top F}{T}\right) \Lambda_t^\top + \hat{S}_u.

5. **PSD Repair:**A final positive semidefinite repair is performed by flooring eigenvalues atfloor_value and symmetrizing the matrix.

Value

A list containing:

• best_rho: The selected shrinkage parameter\hat{\rho}_t for the local residual covariance.

• residual_cov: The shrunk residual covariance matrix\hat{\Sigma}_e(T).

• total_cov: The final estimated time-varying covariance matrix\Sigma_R(t).

• loadings: The local factor loadings\Lambda_t from the local PCA.

• naive_resid_cov: The raw (unshrunk) residual covariance matrix.

#' Expanding Window Time-Varying Minimum Variance Portfolio Optimization

Description

This function performs time-varying minimum variance portfolio (TV-MVP) optimization usingtime-varying covariance estimation based on Local Principal Component Analysis (Local PCA). Theoptimization is performed over a Expanding window, with periodic rebalancing.The procedure is available either as a stand-alone function or as a method inthe 'TVMVP' R6 class.

Usage

expanding_tvmvp(  obj,  initial_window,  rebal_period,  max_factors,  return_type = "daily",  kernel_func = epanechnikov_kernel,  rf = NULL,  M0 = 10,  rho_grid = seq(0.005, 2, length.out = 30),  floor_value = 1e-12,  epsilon2 = 1e-06)

Arguments

Details

Two usage styles:#'

# Function interfaceresults <- expanding_tvmvp(  obj = tv,  initial_window = 50,  rebal_period = 20,  max_factors = 3,  return_type = "daily",  rf = NULL)# R6 method interfacetv <- TVMVP$new()tv$set_data(returns)results <- tv$expanding_tvmvp(  initial_window = 50,  rebal_period = 20,  max_factors = 3,  return_type = "daily",  rf = NULL)

The function implements a Expanding time-varying PCA approach to estimate latent factor structuresand uses a sparse residual covariance estimation method to improve covariance matrix estimation.The covariance matrix is used to determine the global minimum variance portfolio (MVP), which isrebalanced periodically according to the specified 'rebal_period'. The number of factors isdetermined by a BIC-type information criterion using the function 'determine_factors', updatedyearly. The bandwidth is determined by Silverman's rule of thumb, updated each rebalancing period.

If 'rf' is 'NULL', the risk-free rate is assumed to be zero.

Value

An R6 object of classExpandingWindow with the following accessible elements:

summary

A data frame of summary statistics for the TV-MVP and equal-weight portfolios, including cumulative excess return (CER), mean excess returns (MER), Sharpe ratio (SR), and standard deviation (SD) (raw and annualized).

TVMVP

A list containing rebalancing dates, estimated portfolio weights, and excess returns for the TV-MVP strategy.

Equal

A list with similar structure for the equal-weight portfolio.

References

Lillrank, E. (2025).A Time-Varying Factor Approach to Covariance Estimation

Fan, Q., Wu, R., Yang, Y., & Zhong, W. (2024). Time-varying minimum variance portfolio. Journal of Econometrics, 239(2), 105339.

Examples

# Generate random returns for 20 assets over 100 periodsset.seed(123)returns <- matrix(rnorm(20*100), nrow = 100, ncol = 20)# Initialize objecttv <- TVMVP$new()tv$set_data(returns)# Run Expanding TV-MVP optimizationresults <- expanding_tvmvp(  obj = tv,  initial_window = 50,  rebal_period = 20,  max_factors = 3,  return_type = "daily",  kernel_func = epanechnikov_kernel,  rf = NULL)# R6 method interfaceresults <- tv$expanding_tvmvp(  initial_window = 50,  rebal_period = 20,  max_factors = 3,  return_type = "daily",  rf = NULL)# Print summaryprint(results)# Plot cumulative log excess returnsplot(results)

the function will return the size of objand it is smart in the sense that it will choose the suitable unit

Description

the function will return the size of objand it is smart in the sense that it will choose the suitable unit

Usage

get_object_size(obj)

Arguments

Test for Time-Varying Covariance via Local PCA and Bootstrap

Description

This function performs a hypothesis test for time-varying covariance in asset returns based on Su and Wang (2017).It first standardizes the input returns and then computes a time-varying covariance estimatorusing a local principal component analysis (Local PCA) approach. The test statisticJ_{pT}is computed and its significance is assessed using a bootstrap procedure. The procedure is available either as a stand-alonefunction or as a method in the 'TVMVP' R6 class.

Usage

hyptest(returns, m, B = 200, kernel_func = epanechnikov_kernel)

Arguments

Details

Two usage styles:

# Function interfacehyptest(returns, m=2)# R6 method interfacetv <- TVMVP$new()tv$set_data(returns)tv$determine_factors(max_m=5)tv$hyptest()tvtv$get_bootstrap()         # prints bootstrap test statistics

When using the method form, if 'm' are omitted,they default to values stored in the object. Results are cached andretrievable via class methods.

The function follows the steps below:

1. Standardizes the returns.

2. Computes the optimal bandwidth using the Silverman rule.

3. Performs a local PCA on the standardized returns to extract local factors and loadings.

4. Computes a global factor model via singular value decomposition (SVD) to obtain global factors.

5. Calculates residuals by comparing the local PCA reconstruction to the standardized returns.

6. Computes a test statisticJ_{pT} based on a function of the residuals and covariance estimates as:

   \hat{J}_{pT} = \frac{T p^{1/2} h^{1/2} \hat{M} - \hat{\mathbb{B}}_{pT}}{\sqrt{\hat{\mathbb{V}}_{pT}}},

   where:

   \hat{M} = \frac{1}{pT} \sum_{i=1}^p \sum_{t=1}^T \left(\hat{\lambda}_{it}' \hat{F}_t - \tilde{\lambda}_{i0}' \tilde{F}_t\right),

   \hat{\mathbb{B}}_{pT} = \frac{h^{1/2}}{T^2 p^{1/2}} \sum_{i=1}^p \sum_{t=1}^T \sum_{s=1}^T \left(k_{h,st} \hat{F}_s' \hat{F}_t - \tilde{F}_s' \tilde{F}_t\right)^2 \hat{e}_{is}^2,

   and

   \hat{\mathbb{V}}_{pT} = \frac{2}{p h T^2} \sum_{1\leq s \neq r \leq T} \bar{k}_{sr}^2 \left(\hat{F}_s' \hat{\Sigma}_F \hat{F}_r \right)^2 \left(\hat{e}_r' \hat{e}_s \right)^2.

7. A bootstrap procedure is then used to compute the distribution ofJ_{pT} and derive a p-value.

The function prints a message indicating the strength of evidence for time-varying covariance based on the p-value.

Value

A list containing:

References

Su, L., & Wang, X. (2017). On time-varying factor models: Estimation and testing. Journal of Econometrics, 198(1), 84–101

Examples

# Simulate some random returns (e.g., 100 periods, 30 assets)set.seed(123)returns <- matrix(rnorm(100*30, mean = 0, sd = 0.02), nrow = 100, ncol = 30)# Test for time-varying covariance using 3 factors and 10 bootstrap replicationstest_result <- hyptest(returns, m = 3, B = 10, kernel_func = epanechnikov_kernel)# Print test statistic and p-valueprint(test_result$J_NT)print(test_result$p_value)# Or use R6 method interfacetv <- TVMVP$new()tv$set_data(returns)tv$determine_factors(max_m=5)tv$hyptest(iB = 10, kernel_func = epanechnikov_kernel)tvtv$get_bootstrap()         # prints bootstrap test statistics

Perform Local PCA Over Time

Description

This function performs a local principal component analysis (PCA) on asset returns for each time period, aggregating the results over time. It calls an internal functionlocal_pca() at each time index to extract local factors, loadings, and one-step-ahead factor estimates, and stores these results in lists. It uses previously computed factors to align the sign of the new factors.

Usage

localPCA(returns, bandwidth, m, kernel_func = epanechnikov_kernel)

Arguments

Details

The function processes the input returns overT time periods by iteratively calling thelocal_pca() function. For each timet_i:

1. Kernel weights are computed using the specifiedkernel_func andbandwidth.

2. The returns are weighted by the square root of these kernel weights.

3. An eigen decomposition is performed on the weighted returns' covariance matrix to extract the topm eigenvectors, which are scaled by sqrt(T) to form the local factors.

4. The signs of the new factors are aligned with those of the previous factors.

5. The factor loadings are computed by projecting the weighted returns onto the local factors, normalized byT.

6. A second pass computes a one-step-ahead factor estimate for the current time period.

Value

A list with the following components:

• factors: A list of lengthT, where each element is aT × m matrix of local factors.

• loadings: A list of lengthT, where each element is ap × m matrix of factor loadings.

• m: The number of factors extracted.

• weights: A list of lengthT, where each element is a vector of kernel weights used at that time point.

• f_hat: AT × m matrix of one-step-ahead factor estimates.

References

Su, L., & Wang, X. (2017). On time-varying factor models: Estimation and testing. Journal of Econometrics, 198(1), 84–101.

Perform Local Principal Component Analysis

Description

This function performs a local principal component analysis (PCA) on asset returns,weighted by a specified kernel function. It extracts local factors and loadings from the weightedreturns and computes a factor estimate. Optionally, previously estimated factors canbe provided to align the new factors' directions.

Usage

local_pca(returns, r, bandwidth, m, kernel_func, prev_F = NULL)

Arguments

Details

The function operates in the following steps:

1. **Kernel Weight Computation:** For each time pointt = 1, \dots, T, the kernel weight is computed usingboundary_kernel(r, t, T, bandwidth, kernel_func). The weighted returns are given by

   X_r = \text{returns} \circ \sqrt{k_h},

   where\circ denotes element-wise multiplication andk_h is the vector of kernel weights.

2. **Eigen Decomposition:** The function computes the eigen decomposition of the matrixX_r X_r^\top and orders the eigenvalues indescending order. The topm eigenvectors are scaled by\sqrt{T} to form the local factors:

   \hat{F}_r = \sqrt{T} \, \text{eigvecs}_{1:m}.

3. **Direction Alignment:** If previous factors (prev_F) are provided, the function aligns the signs of the new factors with the previous ones by checking the correlation and flipping the sign if the correlation is negative.

4. **Loadings Computation:** The loadings are computed by projecting the weighted returns onto the factors:

   \Lambda_r = \frac{1}{T} X_r^\top \hat{F}_r,

   where the result is transposed to yield ap × m matrix.

5. **One-Step-Ahead Factor Estimation:** A second pass computes the factor estimate for the current time indexr by solving

   \hat{F}_r = \left(\Lambda_r^\top \Lambda_r\right)^{-1} \Lambda_r^\top R_r,

   whereR_r is the return vector at timer.

Value

A list with the following components:

• factors: AT × m matrix of local factors estimated from the weighted returns.

• f_hat: A1 × m vector containing the factor estimate for timer.

• loadings: Ap × m matrix of factor loadings.

• w_r: A numeric vector of kernel weights used in the computation.

Predict Optimal Portfolio Weights Using Time-Varying Covariance Estimation

Description

This function estimates optimal portfolio weights using a time-varying covariance matrixderived from Local Principal Component Analysis (Local PCA). The procedure is available either as a stand-alonefunction or as a method in the 'TVMVP' R6 class. It computes the following portfolios:

1. Global Minimum Variance Portfolio (MVP)

2. Maximum Sharpe Ratio Portfolio (ifmax_SR = TRUE)

3. Return-Constrained Minimum Variance Portfolio (ifmin_return is provided)

Usage

predict_portfolio(  obj,  horizon = 1,  max_factors = 3,  kernel_func = epanechnikov_kernel,  min_return = NULL,  max_SR = NULL,  rf = NULL)

Arguments

Details

Two usage styles:

#'

# R6 method interfacetv <- TVMVP$new()tv$set_data(returns)tv$determine_factors(max_m=5)prediction <- tv$predict_portfolio(horizon = 1, min_return = 0.01, max_SR = TRUE)#' # Function interfaceprediction <- predict_portfolio(obj, horizon = 5, m = 2, min_return = 0.01, max_SR=TRUE)

The methods can then be used onprediction to retrieve the weights.

The function estimates a time-varying covariance matrix using Local PCA:

\hat{\Sigma}_{r,t}=\hat{\Lambda}_t \hat{\Sigma}_F \hat{\Lambda}_t' + \tilde{\Sigma}_e

Where\hat{\Lambda}_t is the factor loadings at time t,\hat{\Sigma}_F is the factor covariance matrix, and\tilde{\Sigma}_e is regularized covariance matrix of the idiosyncratic errors.

It forecasts asset-level expected returns using a simple ARIMA model selection procedure and uses these in portfolio optimization.Optimization strategies include:

• Global minimum variance (analytical)

• Maximum Sharpe ratio (ifmax_SR = TRUE)

• Minimum variance with expected return constraint (ifmin_return is provided)

Value

An object of classPortfolioPredictions (an R6 object) with:

summary

A data frame of evaluation metrics (expected return, risk, Sharpe ratio) for all computed portfolios.

MVP

A list containing the weights, expected return, risk, and Sharpe ratio for the Global Minimum Variance Portfolio.

max_SR

(Optional) A list with metrics for the Maximum Sharpe Ratio Portfolio.

MVPConstrained

(Optional) A list with metrics for the Return-Constrained Portfolio.

Methods

The returned object includes:

• $print(): Nicely prints summary and portfolio access information.

• $getWeights(method = c("MVP", "max_SR", "MVPConstrained")): Retrieves the weights for the selected portfolio.

References

Lillrank, E. (2025).A Time-Varying Factor Approach to Covariance Estimation

Fan, Q., Wu, R., Yang, Y., & Zhong, W. (2024). Time-varying minimum variance portfolio. Journal of Econometrics, 239(2), 105339.

Examples

set.seed(123)returns <- matrix(rnorm(200 * 20, mean = 0, sd = 0.02), ncol = 20)# Initialize objecttv <- TVMVP$new()tv$set_data(returns)# Optimize weights and predict returnsresult <- predict_portfolio(  tv,  horizon = 5,  m = 3,  min_return = 0.02,  max_SR = TRUE)# Print the portfolio performance summaryprint(result)# Access MVP weightsresult$getWeights("MVP")# Access Max Sharpe weights (if computed)result$getWeights("max_SR")# Or use R6 method interfacetv$determine_factors(max_m=5)prediction <- tv$predict_portfolio(horizon = 1, min_return)predictionprediction$getWeights("MVPConstrained")

Estimate Residuals from Factor Model

Description

This function estimates the residuals of asset returns after removing the effectof factor-driven returns.

Usage

residuals(factors, loadings_list, returns)

Arguments

Details

For each time periodt, the function models the asset returns as:

R_t = F_t \Lambda_t + \epsilon_t

whereR_t is the vector of asset returns,F_t is the t'th row of the factor matrix,\Lambda_t is the loadings matrix, and\epsilon_t represents the residuals.

The residuals are computed as the difference between actual returns and the modeledreturns.

Value

A matrix of residuals where each row corresponds to a time period and eachcolumn corresponds to an asset.

Compute Bandwidth Parameter Using Silverman's Rule of Thumb

Description

This function calculates the bandwidth parameter for kernel functions using Silverman's rule of thumb,which is commonly used in kernel density estimation to determine an appropriate bandwidth. The procedure is available either as a stand-alone function or as a method in the 'TVMVP' R6 class.

Usage

silverman(returns)

Arguments

Details

Silverman's rule of thumb for bandwidth selection is given by:

bandwidth = \frac{2.35}{\sqrt{12}} \times T^{-0.2} \times p^{-0.1}

whereT is the number of time periods andp is the number of assets.

Two usage styles:

# Function interfacebw <- silverman(returns)# R6 method interfacetv <- TVMVP$new()tv$set_data(returns)tv$silverman()

Value

A numeric value representing the computed bandwidth parameter based on Silverman's rule.

Examples

# Simulate data for 50 assets over 200 time periodsset.seed(123)T <- 200p <- 50returns <- matrix(rnorm(T * p, mean = 0.001, sd = 0.02), ncol = p)# Compute bandwidth using Silverman's rule of thumbbw <- silverman(returns)print(bw)tv <- TVMVP$new()tv$set_data(returns)tv$silverman()

Compute the Square Root of a Matrix

Description

Computes the square root of a symmetric matrix via eigen decomposition.Negative eigenvalues are handled by taking the square root of their absolute values.

Usage

sqrt_matrix(Amat)

Arguments

Value

A matrix that is the square root ofAmat.

Estimate Time-Varying Covariance Matrix Using Local PCA

Description

This function estimates a time-varying covariance matrix using local principal componentanalysis and the soft thresholding for residual shrinkage. By default, only the totalcovariance matrix is returned. Optionally, the user can retrieve all intermediatecomponents of the estimation process. The procedure is available either as astand-alone function or as a method in the 'TVMVP' R6 class.

Usage

time_varying_cov(  obj,  max_factors = 3,  kernel_func = epanechnikov_kernel,  M0 = 10,  rho_grid = seq(0.005, 2, length.out = 30),  floor_value = 1e-12,  epsilon2 = 1e-06,  full_output = FALSE)

Arguments

Details

The function estimates a time-varying covariance matrix using Local PCA:

\hat{\Sigma}_{r,t}=\hat{\Lambda}_t \hat{\Sigma}_F \hat{\Lambda}_t' + \tilde{\Sigma}_e

Where\hat{\Lambda}_t is the factor loadings at time t,\hat{\Sigma}_F is the factor covariance matrix, and\tilde{\Sigma}_e is regularized covariance matrix of the idiosyncratic errors.

Two usage styles:

# Function interfacetv <- TVMVP$new()tv$set_data(returns)cov <- time_varying_cov(tv, m = 5)# R6 method interfacetv$determine_factors(max_m = 5)cov <- tv$time_varying_cov()

Value

By default, a covariance matrix. Iffull_output = TRUE, a list containing:

• total_cov – the estimated covariance matrix,

• residual_cov – the residual (idiosyncratic) covariance,

• loadings – estimated factor loadings,

• best_rho – optimal shrinkage parameter,

• naive_resid_cov – residual covariance before shrinkage

References

Lillrank, E. (2025).A Time-Varying Factor Approach to Covariance Estimation

Chen, J., Li, D., & Linton, O. (2019). A new semiparametric estimation approach for large dynamic covariance matrices with multiple conditioning variables. Journal of Econometrics, 212(1), 155–176.

Fan, Q., Wu, R., Yang, Y., & Zhong, W. (2024). Time-varying minimum variance portfolio. Journal of Econometrics, 239(2), 105339.

Examples

set.seed(123)returns <- matrix(rnorm(100 * 30), nrow = 100, ncol = 30)# Initialize objecttv <- TVMVP$new()tv$set_data(returns)# Using function interfacetime_varying_cov(obj = tv, m = 3)# Or using R6 methodtv$time_varying_cov(m=3)

Two-Fold Convolution Kernel Function

Description

This function computes the two-fold convolution of a given kernel function with itself.The convolution is evaluated over a range of inputsu and is set to zero outsidethe interval \([-2, 2]\).

Usage

two_fold_convolution_kernel(u, kernel_func)

Arguments

Details

The two-fold convolution kernel is defined as:

K^{(2)}(u) = \int_{-1}^{1} K(v) \cdot K(u - v) \, dv

whereK is the original kernel function. The function evaluates this convolution for eachinputu within the interval \([-2, 2]\) and sets it to zero outside this range.

Value

A numeric vector of two-fold convolution kernel values corresponding to each inputu.