---
title: "IPCA with the Grunfeld dataset"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{IPCA with the Grunfeld dataset}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This vignette walks through `ipca_est()` on the Grunfeld (1958) investment
panel: 11 US firms observed annually from 1935 to 1954. The dataset is
small, balanced, and has only two characteristics, which makes it useful
as a check that the R implementation matches the Python `ipca` package
of Kelly, Pruitt, and Su (2019). The IPCA estimator was developed for
asset returns, but the algorithm only needs a panel and a set of
observable instruments, so the investment-on-fundamentals example
serves as a transparent reference case.

## The IPCA model

Instrumented Principal Components Analysis (IPCA) extracts latent
factors from a panel by letting observable, asset-specific
characteristics act as instruments for otherwise-unobservable
conditional loadings. Without that link the loadings have to be either
pre-specified or estimated as free parameters; IPCA instead ties them
to data the researcher already has.

In the version used here (no intercept, no anomaly term) the model is

$$r_{i,t} = \mathbf{z}_{i,t}^\top \mathbf{\Gamma}_{\!\beta}\, \mathbf{f}_t
          + \varepsilon_{i,t},$$

where $r_{i,t}$ is the outcome of unit $i$ at time $t$ (gross
investment, in the Grunfeld example), $\mathbf{z}_{i,t}$ is the
$L$-vector of observable characteristics, $\mathbf{\Gamma}_{\!\beta}$
is the $L \times K$ matrix that maps characteristics into factor
loadings, and $\mathbf{f}_t$ is the $K$-vector of latent factors. The
unit-specific, time-varying loading is therefore
$\beta_{i,t} = \mathbf{z}_{i,t}^\top \mathbf{\Gamma}_{\!\beta}$, i.e. a
linear projection of characteristics onto the $K$-dimensional factor
space.

Estimation alternates between two ordinary-least-squares steps: solve
for $\mathbf{f}_t$ given $\mathbf{\Gamma}_{\!\beta}$, then update
$\mathbf{\Gamma}_{\!\beta}$ given $\mathbf{f}_t$. After each pass the
solution is rotated by an SVD so that
$\mathbf{\Gamma}_{\!\beta}^\top \mathbf{\Gamma}_{\!\beta} = \mathbf{I}_K$,
which fixes the otherwise-arbitrary scaling and rotation of the
factors.

## Data preparation

The Grunfeld dataset ships with **sdim**. The dependent variable is
gross investment (`invest`); the two characteristics are market value
(`value`) and capital stock (`capital`).

```{r}
library(sdim)

data(grunfeld)
str(grunfeld)
```

`ipca_est()` expects a $T \times N$ outcome matrix and a $T \times N
\times L$ characteristics array. The loop below reshapes the long panel
into that wide form, keeping the year and firm labels on the margins so
the output is easy to read back:

```{r}
firms <- sort(unique(grunfeld$firm))
years <- sort(unique(grunfeld$year))
N  <- length(firms)
TT <- length(years)

ret <- matrix(NA_real_, TT, N, dimnames = list(years, firms))
Z   <- array(NA_real_, dim = c(TT, N, 2),
             dimnames = list(years, firms, c("value", "capital")))

for (i in seq_along(firms)) {

  idx <- grunfeld$firm == firms[i]
  ret[, i]  <- grunfeld$invest[idx]
  Z[, i, 1] <- grunfeld$value[idx]
  Z[, i, 2] <- grunfeld$capital[idx]

}

cat("ret:", nrow(ret), "x", ncol(ret), "\n")
cat("Z:  ", paste(dim(Z), collapse = " x "), "\n")
```

## Fitting IPCA

We start with a single latent factor. With $K = 1$ the loadings
$\beta_{i,t}$ are a single linear combination of `value` and `capital`,
and $f_t$ is one common time-series:

```{r}
fit <- ipca_est(ret, Z, nfac = 1)
print(fit)
summary(fit)
```

The returned object stores the characteristic loadings
($\mathbf{\Gamma}_{\!\beta}$, named `lambda`) and the estimated factor
path:

```{r}
# How each characteristic maps onto the factor
fit$lambda

# Factor realisations over time
data.frame(year = years, factor = fit$factors[, 1])
```

## Validation against the Python `ipca` package

The Python `ipca` package uses the same Grunfeld panel as its built-in
example and runs the identical alternating-least-squares algorithm. With
one factor and no intercept it reports the following loadings and
factor path, which we hard-code below as a reference:

```{r}
py_gamma   <- c(0.99166014, 0.12888046)
py_factors <- c(
  0.1031968381, 0.0884489515, 0.0838496628, 0.0845069923, 0.0722523449,
  0.0995068155, 0.1228840058, 0.1422623752, 0.1197532025, 0.1179724004,
  0.1087561863, 0.1357521189, 0.1579348267, 0.1660545375, 0.1484923276,
  0.1586634303, 0.1596007400, 0.1759379247, 0.1921695585, 0.2111065868
)
```

Factors are identified only up to sign, so before comparing we flip the
R output if the loading vectors are negatively correlated:

```{r}
r_gamma <- as.numeric(fit$lambda)
r_factors <- as.numeric(fit$factors)

if (cor(r_gamma, py_gamma) < 0) {

  r_gamma   <- -r_gamma
  r_factors <- -r_factors

}

cat("Gamma max |diff|:  ", sprintf("%.2e", max(abs(r_gamma - py_gamma))), "\n")
cat("Factor max |diff|: ", sprintf("%.2e", max(abs(r_factors - py_factors))), "\n")
cat("Factor correlation:", sprintf("%.10f", cor(r_factors, py_factors)), "\n")
```

The maximum absolute differences are at numerical-tolerance levels and
the factor correlation is one to ten decimals.

## Multiple factors

The same call extracts more factors by increasing `nfac`. With $K = 2$
each characteristic effectively contributes its own factor dimension
(since there are only two instruments):

```{r}
fit2 <- ipca_est(ret, Z, nfac = 2)
summary(fit2)
```

## References

Grunfeld, Y. (1958). The Determinants of Corporate Investment.
Ph.D. thesis, Department of Economics, University of Chicago.

Kelly, B. T., Pruitt, S., and Su, Y. (2019). Characteristics are
Covariances: A Unified Model of Risk and Return. *Journal of
Financial Economics*, 134(3), 501--524.