## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(sdim)

data(grunfeld)
str(grunfeld)

## -----------------------------------------------------------------------------
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")

## -----------------------------------------------------------------------------
fit <- ipca_est(ret, Z, nfac = 1)
print(fit)
summary(fit)

## -----------------------------------------------------------------------------
# How each characteristic maps onto the factor
fit$lambda

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

## -----------------------------------------------------------------------------
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
)

## -----------------------------------------------------------------------------
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")

## -----------------------------------------------------------------------------
fit2 <- ipca_est(ret, Z, nfac = 2)
summary(fit2)