How to model RVs that depend on time dependent matrix multiplication?

I have X and Y which are both vectors of 30 AR1 random variables with 1 < t < 100 (discrete).

I also have G_t which is a deterministic (sequence already known) 30 x 30 matrix that depends on time t. Lastly, I have Z_t ~ N(X_t - G_t * Y_t , phi). Basically, at time t, each variable in Z_t is distributed with mean as the corresponding entry of X_t minus some linear combination of the Y_ts.

I am unsure how to model this - specifically how do I set this preordered sequence of G_ts and have pymc3 know which G_t to look for and how do I model my observations data if it would be fed in like a time series?

Here’s my code:

    import pymc3 as pm
    import  pandas as pd
   import theano.tensor as T

   observations = pd.read_csv("myObs.csv") 
    with pm.Model() as model:
        beta_week_t = pm.Uniform('beta_week_t', lower = 0 , upper = 100, shape = 30)

        time = 100
        num = 30
        X = pm.AR1('X',beta_week_t,1, shape = (num,time))
        Y = pm.AR1('Y',beta_week_t,1, shape = (num,time))

        Gt = pass ## What do I put here? 

        Z = pm.AR1('outcome', X -,Y), 1, shape = (num, time) , observed = observations )