Best modeling approach for "two-factor varying slope" model?

Question

I'm new to this forum so I hope this question is appropriate. Please let me know if there is anything I can do to improve the question. I simply have a situation in which I am considering the best modeling approach, and would love insight from statisticians.

For a given security numbered $i$, let $r_{i,t}$ be its time-series of its returns. Let $m_t$ be the time-series of the "market's" return over the same time period. Suppose we can categorize all securities in two ways: by their color and by their shape, and that these are not completely independent of each other. Let's say we want to calculate, for each color-shape combination $(c,s)$, the quantity $$\rho\frac{\sigma_{(c,s)}}{\sigma_{m}},$$ where $\rho$ is the correlation between the market return and the securities with color $c$ and shape $s$, and $\sigma_{c,s}$ represents the standard deviation of the $(c,s)$ security returns and $\sigma_{m}$ the standard deviation of the market return.

One simple estimate of this comes from simple linear regression of the $(c,s)$ security returns onto the market returns, since in this regression $$r_{(c,s)} = \beta_{c,m} r_{m} + \epsilon_{c,m}, $$ where $\beta = \rho\frac{\sigma_{(c,s)}}{\sigma_{m}}$. Simple enough!

But suppose we do not have much (or any) data for some combinations $(c,s)$, and yet this is critically important to estimate. A natural way of overcoming this is to assume that each $\beta_{c,s}$ can be factored into a $\beta_{ c}\beta_{s}$ and to leverage the data you have to get good estimates for $\beta_c$ and $\beta_s$ and then "fill in the gap".

Suppose the goal is to find this $\beta_c$ and $\beta_s$. How can this be done?

Approach 1: use a least squares estimator to find the $\beta_c$ and $\beta_s$ that minimize the squared deviations over all observations.

Pro: simple

Con: This seems bias your estimates of $\beta$ toward colors or shapes which covary least with $X$ and $Y$, or which simply have very high variance.

Approach 2: specify a non-nested multilevel model with a two-factor varying slope, where the covariance between $c$ and $s$ - and hence the residual variance for each "bucket", is estimated from the data, but it's not clear to me that this has ever been done or exactly what assumptions on the error terms are made here. My thought is that one might assume that the error terms for each bucket have variances which are normally distributed with variance depending on the covariance of the color and shape, which is a weighted average of the sample covariance and whatever the analogy is to the group-level mean here.

Pro: seems to better capture the dynamics of each bucket, leaning on the information for each color $c$ or $s$ only when the number of samples for that bucket is low.

Con: I don't know if this is possible or has been done.

Question: Am I am approaching this the right way? It seems to me that a multi-level model is exactly the right thing for this situation. I've heard of varying-slope models, but not a "two-factor" varying slope model. And I'm a little unsure of what the precise underlying assumptions of such a model would be. In particular, the way I phrased it begs the question of how to estimate the covariances for shape-color combinations for which there is very little data.

I'm looking for a practical way to overcome this challenge. How might it be done? I would love to do something like this in Stan.

Edit: I edited this question to remove reference to "covariance matrices", which are not really what I am after.

Answer 1

Essentially any model where the kernel of the posterior distribution in log units could be expressed in the Stan language (by incrementing the target keyword), and you can add any prior assumptions you want on the parameters. Getting the Stan MCMC algorithm to sample sufficiently efficiently is another matter, but you have a good chance of that in this situation.

Your Stan program might look something like

data {
  int<lower=1> C;      // number of colors
  int<lower=1> S;      // number of shapes
  int<lower=1> T;      // number of time periods
  vector[T] r_m;       // market returns
  vector[C * S * T] r; // stacked returns
}
parameters {
  vector[C] beta_c;
  vector[S] beta_s;
  real<lower=0> sigma[C,S];
}
model {
  vector[C * S * T] mu; // conditional mean
  vector[C * S * T] sigma_flat; // error sd
  int pos = 1;
  for (c in 1:C) for (s in 1:S) {
    mu[pos:(pos + T - 1)] = beta_c[c] * beta_s[s] * r_m;
    sigma_flat[pos:(pos + T - 1)] = rep_vector(sigma[c,s], T);
    pos = pos + T;
  }
  target += normal_lpdf(r | mu, sigma_flat);  // log-likelihood
  // illustrative priors
  target += normal_lpdf(beta_c | 0, 1);
  target += normal_lpdf(beta_s | 0, 1);
  target += exponential_lpdf(to_vector(sigma) | 1);
}

But you will have to modify it if your panels are unbalanced, if you have better priors, or if you want to relax the assumption that beta_c and beta_s are independent a priori.