A partialCI Guide

Jonas Rende

Christopher Krauss

Matthew Clegg



The partialCI package fits a partial cointegration model to describe a time series. Partial cointegration (PCI) is a weakening of cointegration, allowing for the residual series to contain a mean-reverting and a random walk component. Analytically, this residual series is described by a partially autoregressive process (PAR – see Summers (1986), Poterba and Summers (1988), and Clegg (2015a)), consisting of a stationary AR-process and a random walk. Whereas classic cointegration in the sense of Engle and Granger (1987) requires all shocks to be transient, PCI is more flexible and allows for permanent shocks as well – a realistic assumption across many (macro)economic applications. Even though neither the residual series, nor its mean-reverting and permanent component are directly observable, estimation is still possible in state space – see Brockwell and Davis (2010) and Durbin and Koopman (2012).

The partial cointegration framework

Model definition

Based on Engle and Granger (1987), Clegg and Krauss (2018) define the concept of partial cointegration as follows: : “The components of the vector \(X_t\) are said to be partially cointegrated of order \(d\), \(b\), denoted \(X_t \sim PCI\left(d,b\right)\), if (i) all components of \(X_t\) are \(I\left(d\right)\); (ii) there exists a vector \(\alpha\) so that \(Z_t = \alpha \prime X_t\) and \(Z_t\) can be decomposed as a sum \(Z_t = R_t + M_t\), where \(R_t \sim I\left(d\right)\) and \(M_t \sim I\left(d - b\right)\).”

Let \(Y_t\) denote the target time series and \(X_{j,t}\) the \(j^{th}\) factor time series at time \(t\), where \(j = \lbrace 1, 2, \dots, k \rbrace\). The target time series and the \(k\) factor time series are partially cointegrated, if a parameter vector \(\iota = \left\lbrace\beta_1, \beta_2, \dots, \beta_k, \rho, \sigma_M, \sigma_R, M_0, M_R\right\rbrace\) exists such that the subsequent model equations are satisfied:

\[ \begin{aligned} Y_{t} &= \beta_1 X_{1,t} + \beta_2 X_{2,t} + ... + \beta_k X_{k,t} + W_t \\ W_t &= M_t + R_t \\ M_t &= \rho M_{t-1} + \varepsilon_{M,t}\\ R_t &= R_{t-1} + \varepsilon_{R,t}\\ \varepsilon_{M,t} &\sim \mathcal{N}\left(0, \sigma^2_M\right)\\ \varepsilon_{R,t} &\sim \mathcal{N}\left(0, \sigma^2_R\right)\\ \beta_j \in \mathbb{R}; \rho &\in \left(-1, 1\right);\sigma^2_M, \sigma^2_R \in \mathbb{R}_0^+. \\ \end{aligned} \] Thereby, \(W_t\) denotes the partially autoregressive process, \(R_t\) the permanent component, \(M_t\) the transient component and \(\beta = \lbrace \beta_1, \beta_2, \dots, \beta_k \rbrace\) is the partially cointegrating vector. The permanent component is modeled as a random walk and the transient component as an AR(1)-process with \(AR(1)\)-coefficient \(\rho\). The corresponding error terms \(\varepsilon_{M,t}\) and \(\varepsilon_{R,t}\) are assumed to follow mutually independent, normally distributed white noise processes with mean zero and variances \(\sigma^2_M\) and \(\sigma^2_R\). A key advantage of modeling the cointegrating process as a partially autoregressive process is that we are able to calculate the proportion of variance attributable to mean-reversion (PVMR), defined as (Clegg and Krauss (2018)), \[ R^2_{MR} = \frac{VAR\left[\left(1-B\right)M_t\right]}{VAR\left[\left(1-B\right)W_t\right]} = \frac{2\sigma^2_M}{2\sigma^2_M + \left(1+\rho\right)\sigma^2_R} , \hspace{0.2cm} R^2_{MR} \in \left[0,1\right], \] where \(B\) denotes the backshift operator. The statistic \(R^2_{MR}\) is useful to assess how close the cointegration process is to either a pure random walk \(\left(R^2_{MR} = 0\right)\) or a pure AR(1)-process \(\left(R^2_{MR} = 1\right)\).

State space represenation

The applied state space transformation is in line with Clegg and Krauss (2018). Given that the PAR process \(W_t\) is not observable, we convert the PCI model into the following state space model, consisting of an observation and a state equation:
\[ \begin{align} X_t &= H Z_t \\ Z_t &= FZ_{t-1} + W_t. \end{align} \] Thereby, \(Z_t\) denotes the state which is assumed to be influenced linearly by the state in the last period and a noise term \(W_t\). The matrix \(F\) is assumed to be time invariant. The observable part is denoted by \(X_t\). By assumption, there is a linear dependence between \(X_t\) and \(Z_t\), captured in the time invariant matrix \(H\).

Estimation of a partial cointegration model

Parameters are estimated via the maximum likelihood (ML) method. Using a quasi-Newton algorithm, the ML method searches for the parameters \(\rho\), \(\sigma^2_M\), \(\sigma^2_R\) and the parameter vector \(\beta\) which maximizes the likelihood function of the associated Kalman filter.

A likelihood ratio test routine for partial cointegration}

A time series is classiefied as PCI if we can individually reject both of the subsequent null hypothesis: (i) pure random walk \(\left(\mathcal{H}_0^{R}\right)\) and (ii) pure stationary AR(1)-process \(\left(\mathcal{H}_0^{M}\right)\) – reflecting the cointegration case. The alternative is in both cases unrestricted PCI \(\left(\mathcal{H}_1^{MR}\right)\). Since the proposed null models are all nested in the PCI model (i) and (ii) can be tested with a likelihood ratio test scheme (Neyman and Pearson (1933)). Given that under the null hypothesis we are solely testing parameters at the boundaries of the corresponding parameter space the regularity condition underlying the theorem of Wilks (Wilks (1938)) that the true parameter has to be an inner point of the parameter space does not hold. At boundaries standard asymptotics do not hold. There is a stream of the literature adressing the problem of likelihood ratio tests under such non-standard conditions (Shapiro (1988)). The asymptotic distribution of the likelihood test statistic is a mixture of \(\chi^2\) distributions (Shapiro (1988)): \[\begin{equation}\label{eq: chibar} \sum\limits_{m=0}^k w_m \chi^2_{(m)}, \end{equation}\]

where the weights \(w_i\), with \(i \in \{0, ... , k\}\) sum up to one. Let \(m\) denote the degrees of freedom assosciated with the \(\chi^2\) distribution of interest and \(k\) is equal to the parameter difference of the full and the null model. Note that \(\chi^2_{(0)} \equiv 0\), indicating a point mass at zero. Under Wilks theorem the weight for the \(\chi^2_{(k)}\) distribution is set to one. Therefore, the resulting p-values assuming Wilks theorem holds are more conservative compared to the p-values associated with a mixed distribution. For large sample sizes we calculate p-values under Wilks theorem which can be seen as a conservative approximation of the p-values calculated from the mixed \(\chi^2\) distribution (Stoel et al. (2006)). To appropriately handle small or medium sample sizes p-values are parameterically bootstrapped as proposed by MacKinnon (2009).

To address the multiple testing problem we provide a Bonferroni (Dunn (1961)) and a Holm correction (Holm (1979)).

Using the PCI package

The main functions of the partialCI package are fit.pci(), test.pci(), statehistory.pci(), and hedge.pci().


The function fit.pci() fits a partial cointegration model to a given collection of time series.

fit.pci(Y, X, pci_opt_method = c("jp", "twostep"), par_model = c("par", "ar1", "rw"), robust = FALSE, nu = 5} 


The test.pci() function tests the goodness of fit of a PCI model.

test.pci(Y, X, pci\_opt\_method = c("jp", "twostep"), irobust = FALSE, inu = 5,  null\_hyp = c("par", "rw",  "ar1"), imethod = c("wilk", "boot"), inrep = 999, istart.seed = 1, alpha = 0.05, use.multicore = FALSE)


To estimate the sequence of hidden states the statehistory.pci() function can be applied.

statehistory.pci(A, data = A\$data, basis = A\$basis)}


The function hedge.pci() finds those k factors from a predefined set of factors which yield the best fit to the target time series.

hedge.pci(Y, X, maxfact = 10, use.multicore = TRUE, minimum.stepsize = 0, verbose = TRUE, exclude.cols = c(), search_type = c("lasso", "full", "limited"), pci_opt_method=c("jp", "twostep"))}


As an introductory example, we explore the relationship between Royal Dutch Shell plc A (RDS-A) and Royal Dutch Shell plc B (RDS-B), using daily (closing) price data from 1 January 2006 to 1 January 2018. RDS-A (Royal Dutch Shell plc - A (2018)) and RDS-B (Royal Dutch Shell plc - B (2018)) data are downloaded from Yahoo Finance. To download the price data we use the getYahooData() function, implemented in the R package TTR (Ulrich (2016)).

    RDSA<-getYahooData("RDS-A", 20060101, 20180101)$Close
    RDSB<-getYahooData("RDS-B", 20060101, 20180101)$Close

A classic cointegration analysis yields that the two time series are not cointegrated.


    egcm_finance <- egcm(RDSA,RDSB,include.const = T)

In particular, we apply the two-step approach of Engle and Granger (1987) implemented in the R package egcm (Clegg (2015c)).

The following residual plot (code: plot(egcm\_finance\$residuals,type = "l")) suggests that the residual series is not purely mean-reverting, but rather shows a stochastical trend as well as a mean-reverting behavior.

Residual plot classic cointegration: RDS-A and RDS-B (1.01.2006 - 1.01.2018, daily)

Residual plot classic cointegration: RDS-A and RDS-B (1.01.2006 - 1.01.2018, daily)

Hence, it is not suprising that RDS-A and RDS-B are not cointegrated. Using the PCI framework, we are able to fit a PCI model to RDS-A and RDS-B.

PCI_RDSA_RDSB<-fit.pci(RDSA, RDSB, pci_opt_method = c("jp"), par_model = c("par"), robust = FALSE, nu = 5)

The R output is given as

Fitted values for PCI model
    Y[t] = X[t] %*% beta + M[t] + R[t]
    M[t] = rho * M[t-1] + eps_M [t], eps_M[t] ~ N(0, sigma_M^2)
    R[t] = R[t-1] + eps_R [t], eps_R[t] ~ N(0, sigma_R^2)

            Estimate Std. Err
    beta_Close   0.9232   0.0037
    rho          0.3549   0.0935
    sigma_M      0.0937   0.0070
    sigma_R      0.1085   0.0061

    -LL = -1580.80, R^2[MR] = 0.524,

where beta_Close denotes the partially cointegrating coefficient. The PVMR of 0.52 suggests that the spread time series also exhibits a clear mean-reverting behavior.

In the subsequent step, we utilize the test.pci() function to check whether RDS-A and RDS-B are partially cointegrated. Given the large sample size we will utilize imethod = "wilk".

The R code

test.pci(RDSA, RDSB, alpha = 0.05, null_hyp = "par", imethod = "wilk")

leads to the following output:

  LR test of [RW or CI(1)] vs Almost PCI(1) (joint penalty,wilk)
  data:  RDSA / RDSB

  Hypothesis    Statistic  p-value    alpha alpha_bonf alpha_holm
  Random Walk      -61.55    0.000    0.050    0.025    0.050
  AR(1)            -59.69    0.000    0.050    0.025    0.025.

Recall that a time series is classified as partially cointegrated, if and only if the random walk as well as the AR(1)-hypotheses are rejected. The \(p\)-value of 0.000 for both null model indicates that RDS-A and RDS-B are partially cointegrated in the considered period of time.

Next, we demonstrate the use of the statehistory.pci() function which allows to estimate and extract the hidden states. The R code,


results in the R output:

              Y       Yhat      Z          M         R       eps_M      eps_R
 2006-01-03 32.2123 31.6496 0.562708  0.00000000 0.562708  0.0000000  0.0000000
 2006-01-04 32.5445 31.9223 0.622142  0.01758541 0.604556  0.0175854  0.0418485
 2006-01-05 32.1519 31.6260 0.525826 -0.01890081 0.544726 -0.0251415 -0.0598301
 2006-01-06 32.7659 32.1575 0.608460  0.01413479 0.594325  0.0208423  0.0495991
 2017-12-22 63.0633 59.6017 3.461665 -0.04153269 3.503198 -0.0802318 -0.1909302
 2017-12-26 63.8210 60.1070 3.714061  0.05201243 3.662049  0.0667515  0.1588508
 2017-12-27 63.5237 60.0804 3.443331 -0.05171763 3.495049 -0.0701758 -0.1669996
 2017-12-28 63.7827 60.2666 3.516118 -0.00668903 3.522807  0.0116645  0.0277584
 2017-12-29 63.9841 60.5414 3.442713 -0.02536998 3.468083 -0.0229961 -0.0547247.

The latter table covers the estimates of the hidden states \(M\) and \(R\) as well as the corresponding error terms eps_M and eps_R. Z is equal to the sum of \(M\) and \(R\). The estimate of the target time series is denoted by Yhat.

The subsequent figure illustrates a plot of the extracted mean-reverting component of the spread associated with the RDS-A and RDS-B price time series (`plot(statehistory.pci(PCI_RDSA_RDSB)[,4],type = "l",ylab = "", xlab = "")).

Mean-reverting component RDS-A and RDS-B (1.01.2006 - 1.01.2018, daily)

Mean-reverting component RDS-A and RDS-B (1.01.2006 - 1.01.2018, daily)

The horizontal blue lines are equal to two times the historical standard deviation in absolute terms of the mean-reverting component. A pairs trading strategy could exploit the mean-reverting behavior of \(M_t\). Note that this example is in-sample; for a true out-of-sample application see Clegg and Krauss (2018).

We continue with using hedge.pci() to find the set of sector ETFs forming the best hedging portfolio for the SPY index (S&P500 index).

Thereby, the R code,

    sectorETFS <- c("XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY")
    prices <- multigetYahooPrices(c("SPY", sectorETFS), start=20060101)

    hedge.pci(prices[,"SPY"], prices),

results in the subsequent output:

-LL     LR[rw]  p[rw]   p[mr]  rho  R^2[MR] Fac.| Factor coefficients
2521.2 -19.281  0.000  0.000  0.593  0.407  XLI | 3.033 
1891.8 -41.011  0.000  0.000  0.350  0.435  XLY | 1.848 1.182 
1558.8 -47.619  0.000  0.000  0.350  0.461  XLV | 1.627 0.915 0.679 
1302.6 -52.992  0.000  0.000  0.350  0.483  XLU | 1.490 0.895 0.596 0.583 
 724.5 -64.165  0.000  0.000  0.382  0.552  XLK | 1.209 0.544 0.479 0.491 1.265.  

The table summarizes information about the best hedging portfolio, where each row corresponds to an increasing number of factors. Row 1: The best single-factor hedging portfolio comprises XLI (industrials) as only factor. Row 2: The best two-factor hedging portfolio consists of XLI and XLY (consumer discretionary). As such, XLY leads to the best improvement of the LRT score among all remaining factors. Row 3 includes XLV (health care) for the three-factor portfolio and row 4 XLK (technology) for the best four-factor portfolio. The last row corresponds to the overall best fit out of the nine potential sector ETFs, based on the LRT score. Note that for all rows, the union of random walk and AR(1)-null hypothesis is rejected at the 5 percent significant level, so we find a PCI model at each step.


Brockwell, P. J., and R. A. Davis. 2010. Introduction to time series and forecasting. 2. ed., corr. at 8th print. Springer texts in statistics. New York [u.a.]: Springer.

Clegg, M. 2015a. “Modeling time series with both permanent and transient components using the partially autoregressive model.” SSRN Electronic Journal. http://dx.doi.org/10.2139/ssrn.2556957.

———. 2015c. “egcm: Engle-Granger cointegration models.” https://CRAN.R-project.org/package=egcm.

Clegg, M., and C. Krauss. 2018. “Pairs trading with partial cointegration.” Quantitative Finance 18 (1): 121–38.

Dunn, Olive Jean. 1961. “Multiple Comparisons Among Means.” American Statistical Association 56 (293): 52–64.

Durbin, J., and S. J. Koopman. 2012. Time series analysis by state space methods. Second. Vol. 38. Oxford statistical science series. Oxford: Oxford University Press.

Engle, R. F., and C. W. J. Granger. 1987. “Co-Integration and error correction: Representation, estimation, and testing.” Econometrica 55 (2): 251.

Friedman, J., T. Hastie, and R. Tibshirani. 2010. “Regularization paths for generalized linear models via coordinate descent.” Journal of Statistical Software 33 (1): 1–22. https://CRAN.R-project.org/package=glmnet.

Holm, Sture. 1979. “A Simple Sequentially Rejective Multiple Test Procedure.” Scandinavian Journal of Statistics 6 (2): 65–70.

MacKinnon, James G. 2009. “Bootstrap Hypothesis Testing: 6.” In Handbook of Computational Econometrics, 183–213. Wiley-Blackwell.

Neyman, J., and E. S. Pearson. 1933. “On the Problem of the Most Efficient Tests of Statistical Hypotheses.” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 231: 289–337.

Poterba, J. M., and L. H. Summers. 1988. “Mean reversion in stock prices.” Journal of Financial Economics 22 (1): 27–59.

Royal Dutch Shell plc - A. 2018. Historical data. retrieved from Yahoo!Finance. https://finance.yahoo.com/quote/RDS-A/history?p=RDS-A.

Royal Dutch Shell plc - B. 2018. Historical data. retrieved from Yahoo!Finance. https://finance.yahoo.com/quote/RDS-B/history?p=RDS-B.

Shapiro, A. 1988. “Towards a Unified Theory of Inequality Constrained Testing in Multivariate Analysis.” International Statistical Review / Revue Internationale de Statistique 56 (1): 49.

Stoel, Reinoud D., Francisca Galindo Garre, Conor Dolan, and Godfried van den Wittenboer. 2006. “On the likelihood ratio test in structural equation modeling when parameters are subject to boundary constraints.” Psychological methods 11 (4): 439–55.

Summers, L. H. 1986. “Does the stock market rationally reflect fundamental values?” The Journal of Finance 41 (3): 591.

Ulrich, J. 2016. “TTR: Technical Trading Rules.” https://CRAN.R-project.org/package=TTR.

Wilks, S. S. 1938. “The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses.” Ann. Math. Statist. 9 (1): 60–62.