From a4eb9a91741d4a83c6fb1f6b20dafee6a74ed316 Mon Sep 17 00:00:00 2001
From: Christoph Rust <christoph.g.rust@gmail.com>
Date: Thu, 13 Feb 2025 10:54:37 +0100
Subject: [PATCH] feat: change dpca's default values

Matteo Barigozzi suggested to use different default values for
bandwidth and evalaution points in the frequency domain of the
dynamic spectral decomposition. This commit applies those default
values.
---
 DESCRIPTION |  2 +-
 R/dpca.R    | 20 ++++++++++++--------
 man/dpca.Rd | 14 +++++++-------
 3 files changed, 20 insertions(+), 16 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 4cb21d8..28c6f0f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,7 +25,7 @@ License: GPL-3 + file LICENSE
 URL: https://github.com/christophrust/dpca
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Suggests:
     testthat (>= 3.0.0),
     freqdom (>= 2.0.3),
diff --git a/R/dpca.R b/R/dpca.R
index 9e68c9e..57ed4f8 100644
--- a/R/dpca.R
+++ b/R/dpca.R
@@ -20,8 +20,8 @@
 #' spectral density is evaluated.
 #'
 #' @param bandwidth Single integer, giving the width of the
-#' lag window estimator. If unspecified, the cube root of the number of time
-#' observations is used as default.
+#' lag window estimator. If unspecified, 0.75 times the square root of the number
+#' of time observations is used as default.
 #'
 #' @param weights Kernel used for the lag window estimation of spectrum.
 #'
@@ -30,6 +30,7 @@
 #'
 #' @param q_max Maximum numer of dynamic factors considered in the data-driven
 #' selection.
+#'
 #' @param n_path Integer vector specifying which (nested) subsets of the
 #' cross section are used in the Hallin & Liska procedure. If unspecified,
 #' a regular sequence of length \code{20} from \code{n/2} to \code{n} is used.
@@ -110,23 +111,24 @@
 dpca <- function(
     x,
     q,
-    freqs = -20:20 / 20 * pi,
-    bandwidth = floor(ncol(x)^(1 / 3)),
+    freqs = -bandwidth:bandwidth / bandwidth * pi,
+    bandwidth = floor(0.75 * sqrt(tx)),
     weights = c(
       "bartlett", "trunc", "tukey", "parzen",
       "bohman", "daniell", "parzen_cogburn_davis"
     ),
     qsel_crit = c("IC1", "IC2"),
     q_max = 15,
-    n_path = floor(seq(nrow(x) / 2, nrow(x), nrow(x) / 20)),
+    n_path = floor(seq(nx / 2, nx, nx / 20)),
     penalties = (
-      bandwidth^(-2) + sqrt(bandwidth / ncol(x)) + 1 / n_path
-    ) * log(pmin(n_path, bandwidth^2, sqrt(ncol(x) / bandwidth))),
+      bandwidth^(-2) + sqrt(bandwidth / tx) + 1 / n_path
+    ) * log(pmin(n_path, bandwidth^2, sqrt(tx / bandwidth))),
     penalty_scales = seq(0, 2, by = 0.01)) {
   if (length(weights) > 1) {
     weights <- "bartlett"
   }
 
+  ## our internal object is n times t
   x <- if (is.ts(x) || "zoo" %in% class(x)) {
     t(x)
   } else if (is.matrix(x)) {
@@ -134,6 +136,8 @@ dpca <- function(
   } else {
     stop("x must either a \"ts\" or \"zoo\" object or a matrix!")
   }
+  nx <- nrow(x)
+  tx <- ncol(x)
 
   if (!missing(q) && (length(q) > 1 || floor(abs(q)) != q)) {
     stop("\"q\" has to be a single positive integer!")
@@ -186,7 +190,7 @@ dpca <- function(
     as.integer(q_max),
     select_q,
     as.integer(n_path),
-    as.integer(rep(ncol(x), length(n_path))),
+    as.integer(rep(tx, length(n_path))),
     as.numeric(penalties),
     as.numeric(penalty_scales),
     PACKAGE = "dpca"
diff --git a/man/dpca.Rd b/man/dpca.Rd
index 6c73944..af49575 100644
--- a/man/dpca.Rd
+++ b/man/dpca.Rd
@@ -7,15 +7,15 @@
 dpca(
   x,
   q,
-  freqs = -20:20/20 * pi,
-  bandwidth = floor(ncol(x)^(1/3)),
+  freqs = -bandwidth:bandwidth/bandwidth * pi,
+  bandwidth = floor(0.75 * sqrt(tx)),
   weights = c("bartlett", "trunc", "tukey", "parzen", "bohman", "daniell",
     "parzen_cogburn_davis"),
   qsel_crit = c("IC1", "IC2"),
   q_max = 15,
-  n_path = floor(seq(nrow(x)/2, nrow(x), nrow(x)/20)),
-  penalties = (bandwidth^(-2) + sqrt(bandwidth/ncol(x)) + 1/n_path) * log(pmin(n_path,
-    bandwidth^2, sqrt(ncol(x)/bandwidth))),
+  n_path = floor(seq(nx/2, nx, nx/20)),
+  penalties = (bandwidth^(-2) + sqrt(bandwidth/tx) + 1/n_path) * log(pmin(n_path,
+    bandwidth^2, sqrt(tx/bandwidth))),
   penalty_scales = seq(0, 2, by = 0.01)
 )
 }
@@ -31,8 +31,8 @@ it is chosed data-driven by using the criterion of Hallin & Liska (2007).}
 spectral density is evaluated.}
 
 \item{bandwidth}{Single integer, giving the width of the
-lag window estimator. If unspecified, the cube root of the number of time
-observations is used as default.}
+lag window estimator. If unspecified, 0.75 times the square root of the number
+of time observations is used as default.}
 
 \item{weights}{Kernel used for the lag window estimation of spectrum.}