How to Create Multivariate AR(1)

I am attempting to create a time-dependent model which depends on an AR1. The point of the model is to model the number of three-pointers a team makes as a gaussian sampled from an AR1. Ideally I’d like to expand this to all 30 NBA teams. If I model just one team, it appears to work quite well.

with pm.Model() as warr_model:
       warr_tau = pm.HalfNormal('warr_tau',sigma = 10)

       warr_coef = pm.Normal('warr_coef',mu = 1, sigma = 20)

       warr_ar1 = pm.AR1('warr_ar1',k = warr_coef, tau_e = warr_tau, shape = len(warr_fg3m))![21%20pm|690x266](upload://yHYqSwFRqNf4GkfNM5mIHTh3ATa.png) 

       game_sd = pm.HalfNormal("game_sd", sigma = 20)

       y = pm.Normal('y',mu = warr_ar1,sigma = game_sd, observed = warr_fg3m)

I get normal results, as seen right here

If I do someething similar with New Orleans I also get normal results:

However, when I try to combine these two into a multivariate AR1 I get wacky results. Note in the code below, the results of New Orleans and the Warriors is independent.

with pm.Model() as warr_nola:
# This iss a combined model where we are doing analysis on both NOLA and Warriors
        comb_tau = pm.HalfNormal('comb_tau',sigma = 10,shape =2)

       comb_coef = pm.Normal('comb_coef',mu = 1, sigma = 20, shape = 2 )

       comb_ar1 = pm.AR1('comb_ar1',k = comb_coef, tau_e = comb_tau, shape = (40,2))

       comb_sd = pm.HalfNormal("comb_sd", sigma = 20,shape = 2)

       y_comb = pm.Normal('y_comb',mu = comb_ar1,sigma = comb_sd, observed = np.stack([nola_fg3m,warr_fg3m],axis =1))

If I sample from this distribution, I get pretty wacky results that are significantly different:

Can someone please explain what’s going wrong. Also the wierdness is combing from the fact that my k-coefficient goes from ~.98 to ~.02 when I change from a univariate to a multivariate case.