# Advanced Data-driven Nonparametric Regression for the Trend in Equidistant Time Series This function runs an iterative plug-in algorithm to find the optimal bandwidth for the estimation of the nonparametric trend in equidistant time series (with short-memory errors) and then employs the resulting bandwidth via either local polynomial or kernel regression. This function allows for more flexibility in its arguments than **msmooth**. ```r tsmooth( y, p = c(1, 3), mu = c(0, 1, 2, 3), Mcf = c("NP", "ARMA", "AR", "MA"), InfR = c("Opt", "Nai", "Var"), bStart = 0.15, bvc = c("Y", "N"), bb = c(0, 1), cb = 0.05, method = c("lpr", "kr") ) ``` ## Arguments - `y`: a numeric vector that contains the time series ordered from past to present. - `p`: an integer `1` (local linear regression) or `3` (local cubic regression); represents the order of polynomial within the local polynomial regression (see also the 'Details' section); is set to `1` by default; is automatically set to `1` if `method = "kr"`. - `mu`: an integer `0`, ..., `3` that represents the smoothness parameter of the kernel weighting function and thus defines the kernel function that will be used within the local polynomial regression; is set to `1` by default. ||| |:-:|:--| |Number|Kernel| |`0`|Uniform Kernel| |`1`|Epanechnikov Kernel| |`2`|Bisquare Kernel| |`3`|Triweight Kernel| - `Mcf`: method for estimating the variance factor $c_f$ by the IPI (see also the 'Details' section); is set to `NP` by default. ||| |:-:|:--| |Method|Explanation| |`"NP"`|Nonparametric estimation using the Bartlett window| |`"ARMA"`|Estimation on the assumption that the residuals follow an ARMA model| |`"AR"`|Estimation on the assumption that the residuals follow an AR model| |`"MA"`|Estimation on the assumption that the residuals follow an MA model| - `InfR`: a character object that represents the inflation rate in the form $h_d = h^a$ for the bandwidth in the estimation of $I[m^(k)]$ (see also the 'Details' section); is set to `"Opt"` by default. ||| |:-:|:--| |Inflation rate|Description| |`"Opt"`|Optimal inflation rate $a_[p,O]$ ( $5/7$ for `p = 1` ; $9/11$ for `p = 3` )| |`"Nai"`|Naive inflation rate $a_[p,N]$ ( $5/9$ for `p = 1` ; $9/13$ for `p = 3` )| |`"Var"`|Stable inflation rate $a_[p,S]$ ( $1/2$ for `p = 1` and `p = 3` )| - `bStart`: a numeric object that indicates the starting value of the bandwidth for the iterative process; should be $> 0$; is set to `0.15` by default. - `bvc`: a character object that indicates whether an enlarged bandwidth is being used for the estimation of the variance factor $c_f$ (see also the 'Details' section) or not; is set to `"Y"` by default. ||| |:-:|:--| |Bandwidth|Description| |**"Y"**|Using an enlarged bandwidth| |**"N"**|Using a bandwidth without enlargement| - `bb`: can be set to `0` or `1`; the parameter controlling the bandwidth used at the boundary; is set to `1` by default. ||| |:-:|:--| |Number (‘bb’)|Estimation procedure at boundary points| |`0`|Fixed bandwidth on one side with possible large bandwidth on the other side at the boundary| |`1`|The $k$ -nearest neighbor method will be used| - `cb`: a numeric value that indicates the percentage of omitted observations on each side of the observation period for the automated bandwidth selection; is set to `0.05` by default. - `method`: the final smoothing approach; `"lpr"` represents the local polynomial regression, whereas `"kr"` implements a kernel regression; is set to `"lpr"` by default. ## Returns The function returns a list with different components: - **AR.BIC**: the Bayesian Information Criterion of the optimal AR($p$) model when estimating the variance factor via autoregressive models (if calculated; calculated for `alg = "OA"` and `alg = "NA"`). - **ARMA.BIC**: the Bayesian Information Criterion of the optimal ARMA($p,q$) model when estimating the variance factor via autoregressive-moving-average models (if calculated; calculated for `alg = "OAM"` and `alg = "NAM"`). - **cb**: the percentage of omitted observations on each side of the observation period. - **b0**: the optimal bandwidth chosen by the IPI-algorithm. - **bb**: the boundary bandwidth method used within the IPI; always equal to 1. - **bStart**: the starting value of the (relative) bandwidth; input argument. - **bvc**: indicates whether an enlarged bandwidth was used for the variance factor estimation or not; depends on the chosen algorithm. - **cf0**: the estimated variance factor; in contrast to the definitions given in the **Details** section, this object actually contains an estimated value of $2\pi c_f$, i.e. it corresponds to the estimated sum of autocovariances. - **cf0.AR**: the estimated variance factor obtained by estimation of autoregressive models (if calculated; `alg = "OA"` or `"NA"`). - **cf0.ARMA**: the estimated variance factor obtained by estimation of autoregressive-moving-average models (if calculated; calculated for `alg = "OAM"` and `alg = "NAM"`). - **cf0.LW**: the estimated variance factor obtained by Lag-Window Spectral Density Estimation following Bühlmann (1996) (if calculated; calculated for algorithms `"A"`, `"B"`, `"O"` and `"N"`). - **cf0.MA**: the estimated variance factor obtained by estimation of moving-average models (if calculated; calculated for `alg = "OM"` and `alg = "NM"`). - **I2**: the estimated value of $I[m^(k)]$. - **InfR**: the setting for the inflation rate according to the chosen algorithm. - **iterations**: the bandwidths of the single iterations steps - **L0.opt**: the optimal bandwidth for the lag-window spectral density estimation (if calculated; calculated for algorithms `"A"`, `"B"`, `"O"` and `"N"`). - **MA.BIC**: the Bayesian Information Criterion of the optimal MA($q$) model when estimating the variance factor via moving-average models (if calculated; calculated for `alg = "OM"` and `alg = "NM"`). - **Mcf**: the estimation method for the variance factor estimation; depends on the chosen algorithm. - **mu**: the smoothness parameter of the second order kernel; input argument. - **n**: the number of observations. - **niterations**: the total number of iterations until convergence. - **orig**: the original input series; input argument. - **p.BIC**: the order p of the optimal AR($p$) or ARMA($p,q$) model when estimating the variance factor via autoregressive or autoregressive-moving average models (if calculated; calculated for `alg = "OA"`, `alg = "NA"`, `alg = "OAM"` and `alg = "NAM"`). - **p**: the order of polynomial used in the IPI-algorithm; also used for the final smoothing, if `method = "lpr"`; input argument. - **q.BIC**: the order $q$ of the optimal MA($q$) or ARMA($p,q$) model when estimating the variance factor via moving-average or autoregressive-moving average models (if calculated; calculated for `alg = "OM"`, `alg = "NM"`, `alg = "OAM"` and `alg = "NAM"`). - **res**: the estimated residual series. - **v**: the considered order of derivative of the trend; is always zero for this function. - **ws**: the weighting system matrix used within the local polynomial regression; this matrix is a condensed version of a complete weighting system matrix; in each row of `ws`, the weights for conducting the smoothing procedure at a specific observation time point can be found; the first $[nb + 0.5]$ rows, where $n$ corresponds to the number of observations, $b$ is the bandwidth considered for smoothing and $[.]$ denotes the integer part, contain the weights at the $[nb + 0.5]$ left-hand boundary points; the weights in row $[nb + 0.5] + 1$ are representative for the estimation at all interior points and the remaining rows contain the weights for the right-hand boundary points; each row has exactly $2[nb + 0.5] + 1$ elements, more specifically the weights for observations of the nearest $2[nb + 0.5] + 1$ time points; moreover, the weights are normalized, i.e. the weights are obtained under consideration of the time points $x_t = t/n$, where $t = 1, 2, ..., n$. - **ye**: the nonparametric estimates of the trend. ## Details The trend is estimated based on the additive nonparametric regression model for an equidistant time series $$ y_t = m(x_t) + \epsilon_t, $$ where $y_t$ is the observed time series, $x_t$ is the rescaled time on the interval $[0, 1]$, $m(x_t)$ is a smooth and deterministic trend function and $\epsilon_t$ are stationary errors with $E(\epsilon_t) = 0$ and short-range dependence (see also Beran and Feng, 2002). With this function $m(x_t)$ can be estimated without a parametric model assumption for the error series. Thus, after estimating and removing the trend, any suitable parametric model, e.g. an ARMA($p,q$) model, can be fitted to the residuals (see `arima`). The iterative-plug-in (IPI) algorithm, which numerically minimizes the Asymptotic Mean Squared Error (AMISE), was proposed by Feng, Gries and Fritz (2020). Define c("$I[m^(k)] =\n$", "$int_[c_b]^[d_b] {m^(k)}^2 dx$"), $\beta_(\nu, k) = int_[-1]^[1] u^k K_(\nu, k)(u) du$ and c("$R(K) = int_[-1]^[1]\n$", "${K_(\\nu, k)}^2 du$"), where $p$ is the order of the polynomial, $k = p + 1$ is the order of the asymptotically equivalent kernel, $\nu$ is the order of the trend function's derivative, $0 \le c_b < d_b \le 1$, $c_f$ is the variance factor and $K_(\nu, k)(u)$ the $k$-th order equivalent kernel obtained for the estimation of $m^(\nu)$ in the interior. $m^(\nu)$ is the $\nu$-th order derivative (c("$\\nu = 0,\n$", "$1, 2, ...$")) of the nonparametric trend. Furthermore, we define $$ C_{1} = \frac{I[m^{(k)}] \beta_{(\nu, k)}^2}{(k!)^2}C_1 =[I[m^(k)]{\beta_(\nu, k)}^2] / (k!)^2 $$ and $$ C_{2} = \frac{2 \pi c_{f} (d_b - c_b) R(K)}{nh^{2 \nu + 1}}C_2 =2\pi(d_b - c_b)R(K)c_f / (nh^[2\nu + 1]) $$ with $h$ being the bandwidth and $n$ being the number of observations. The AMISE is then $$ AMISE(h) = h^{2(k-\nu)}C_{1} + C_{2}.AMISE(h) = h^[2(k - \nu)] C_1 +C_2. $$ The function calculates suitable estimates for $c_f$, the variance factor, and $I[m^(k)]$ over different iterations. In each iteration, a bandwidth is obtained in accordance with the AMISE that once more serves as an input for the following iteration. The process repeats until either convergence or the 40th iteration is reached. For further details on the asymptotic theory or the algorithm, please consult Feng, Gries and Fritz (2020) or Feng et al. (2019). To apply the function, more arguments are needed compared to the similar function `msmooth`: a data input `y`, an order of polynomial `p`, a kernel weighting function defined by the smoothness parameter `mu`, a variance factor estimation method `Mcf`, an inflation rate setting `InfR` (see also Beran and Feng, 2002), a starting value for the relative bandwidth `bStart`, an inflation setting for the variance factor estimation `bvc`, a boundary method `bb`, a boundary cut-off percentage `cb` and a final smoothing method `method`. In fact, aside from the input vector `y`, every argument has a default setting that can be adjusted for the individual case. Theoretically, the initial bandwidth does not affect the selected optimal bandwidth. However, in practice local minima of the AMISE might exist and influence the selected bandwidth. Therefore, the default setting is `bStart = 0.15`, which is a compromise between the starting values `bStart = 0.1` for `p = 1` and `bStart = 0.2` for `p = 3` that were proposed by Feng, Gries and Fritz (2020). In the rare case of a clearly unsuitable optimal bandwidth, a starting bandwidth that differs from the default value is a first possible approach to obtain a better result. Other argument adjustments can be tried as well. For more specific information on the input arguments consult the section **Arguments**. When applying the function, an optimal bandwidth is obtained based on the IPI algorithm proposed by Feng, Gries and Fritz (2020). In a second step, the nonparametric trend of the series is calculated with respect to the chosen bandwidth and the selected regression method (`lpf` or `kr`). Please note that `method = "lpf"` is strongly recommended by the authors. Moreover, it is notable that `p` is automatically set to `1` for `method = "kr"`. The output object is then a list that contains, among other components, the original time series, the estimated trend values and the series without the trend. The default print method for this function delivers only key numbers such as the iteration steps and the generated optimal bandwidth rounded to the fourth decimal. The exact numbers and results such as the estimated nonparametric trend series are saved within the output object and can be addressed via the `$` sign. NOTE: With package version 1.1.0, this function implements C++ code by means of the `Rcpp` and `RcppArmadillo` packages for better performance. ## Examples ```r ### Example 1: US-GDP ### # Logarithm of test data # -> the logarithm of the data is assumed to follow the additive model test_data <- gdpUS y <- log(test_data$GDP) # Applied tsmooth function for the trend result <- tsmooth(y, p = 1, mu = 1, Mcf = "NP", InfR = "Opt", bStart = 0.1, bvc = "Y") trend1 <- result$ye # Plot of the results t <- seq(from = 1947, to = 2019.25, by = 0.25) plot(t, y, type = "l", xlab = "Year", ylab = "log(US-GDP)", bty = "n", lwd = 1, lty = 3, main = "Estimated trend for log-quarterly US-GDP, Q1 1947 - Q2 2019") points(t, trend1, type = "l", col = "red", lwd = 1) title(sub = expression(italic("Figure 1")), col.sub = "gray47", cex.sub = 0.6, adj = 0) result ## Not run: ### Example 2: German Stock Index ### # The following procedure can be considered, if (log-)returns are assumed # to follow a model from the general class of semiparametric GARCH-type # models (including Semi-GARCH, Semi-Log-GARCH and Semi-APARCH models among # others) with a slowly changing variance over time due to a deterministic, # nonparametric scale function. # Obtain the logarithm of the squared returns returns <- diff(log(dax$Close)) # (log-)returns rt <- returns - mean(returns) # demeaned (log-)returns yt <- log(rt^2) # logarithm of the squared returns # Apply 'smoots' function to the log-data, because the logarithm of # the squared returns follows an additive model with a nonparametric trend # function, if the returns are assumed to follow a semiparametric GARCH-type # model. # In this case, the optimal inflation rate is used for p = 3. est <- tsmooth(yt, p = 3, InfR = "Opt") m_xt <- est$ye # estimated trend values # Obtain the standardized returns 'eps' and the scale function 'scale.f' res <- est$res # the detrended log-data C <- -log(mean(exp(res))) # an estimate of a constant value needed # for the retransformation scale.f <- exp((m_xt - C) / 2) # estimated values of the scale function in # the returns eps <- rt / scale.f # the estimated standardized returns # -> 'eps' can now be analyzed by any suitable GARCH-type model. # The total volatilities are then the product of the conditional # volatilities obtained from 'eps' and the scale function 'scale.f'. ## End(Not run) ``` ## References Beran, J. and Feng, Y. (2002). Local polynomial fitting with long-memory, short-memory and antipersistent errors. Annals of the Institute of Statistical Mathematics, 54(2), 291-311. Bühlmann, P. (1996). Locally adaptive lag-window spectral estimation. Journal of Time Series Analysis, 17(3), 247-270. Feng, Y., Gries, T. and Fritz, M. (2020). Data-driven local polynomial for the trend and its derivatives in economic time series. Journal of Nonparametric Statistics, 32:2, 510-533. Feng, Y., Gries, T., Letmathe, S. and Schulz, D. (2019). The smoots package in R for semiparametric modeling of trend stationary time series. Discussion Paper. Paderborn University. Unpublished. ## Author(s) * Yuanhua Feng (Department of Economics, Paderborn University), Author of the Algorithms Website: [https://wiwi.uni-paderborn.de/en/dep4/feng/](https://wiwi.uni-paderborn.de/en/dep4/feng/) * Dominik Schulz (Research Assistant) (Department of Economics, Paderborn University), Package Creator and Maintainer