Forecasting U.S. presidential elections. Example in BDA book

Hi All, I’m just getting started with Stan and I’m having a little trouble. I’m trying to replicate the example forecasting U.S. presidential elections of section 15.2 of Bayesian Data Analysis book. Briefly, the goal of the study is to fit a model for the proportion of Democratic votes in each state from each of the 11 presidential elections from 1948 through 1988.
The model for state s, region r and year t is:

y_{st} \sim N(X_{st}\beta+\gamma_{r(s)t}+\delta_t, \sigma^2) \\ \gamma_{rt} \sim \begin{cases} N(0,\tau_{\gamma1}^2) & \text{for} &r=1,2,3 \\ N(0,\tau_{\gamma2}^2) & \text{for} &r=4 \end{cases} \\ \delta_t \sim N(0, \tau_{\delta}^2)

with a uniform hyperprior distribution on \beta, \sigma, \tau_{\gamma1} , \tau_{\gamma2}, \tau_{\delta}. The study consider 74 predictors: 4 nationwide variables, 9 statewide variables, 6 regional variables, 11 indicator variables for the year and 4 x 11 indicator variables for each combination of year and region.
The data are available here.

I was able to implement only the varying intercept for the years. I can’t figure out to implement the varying intercept for the combination of year and region.

data {
  int<lower=0> N; // number of records
  int<lower=0> K; // number of predictors
  vector[N] y;  //  outcomes (Democratic votes)
  matrix[N, K] x;  // matrix of predictors 
  int<lower=1, upper=11> year[N];  // indicator vector for year

parameters {
  real alpha;
  vector[K] beta;
  vector[11] gamma;
  real<lower=0> tau; 
  real<lower=0> sigma;

model {
  // prior for gamma
  gamma ~ normal(0, tau);
  // likelihood
  y ~ normal(alpha + x*beta + gamma[year], sigma);

I hope the question is appropriate and well-specified for this forum.
Thank you in advance for your help.
Have a good day!!