Discrete Markov Chain and tracking population average

Question

I would appreciate your help.

Suppose a large number of population $N$, where every individual $i$ is attached a fixed number $x_i$. Now every individual switches across two states: blue and red. The switching is a Markov Chain with the transition probability matrix $P$. Now I am interested in tracking the average of $x_i$ across those two groups. Particularly, I am looking for a neat derivation of how the population means vary across time.

Since the population is large I expect that the population weights will evolve as $w_{t-1}'P=w_{t}'$, where $w_t'=[w^{\text{red}} \text{ } w^{\text{blue}}]$. Thus, the probability to be in either state will converge to population weights. But what about the expected value of $X^{\text{red}}_t=\sum x_i$ belonging to the group 1 at time $t$?

My expectation that it can be defined $$w^{\text{red}}_{t} X^{\text{red}}_t = P(\text{stay red}) w^{\text{red}}_{t-1} X^{\text{red}}_{t-1} + P(\text{switch from blue to red}) w^{\text{blue}}_{t-1} X^{\text{blue}}_{t-1}$$

I created a Matlab simulation to illustrate my logic. The simulation agrees with my calculation. But I still would like to rigorously derive the result.

Am I on the right track?

First timer, sorry if something is unclear

Thanks!!!

clear; clc;

% transition_probabilities = [0.1 0.9;0.8 0.2];
transition_probabilities = [0.7 0.2 0.1; 0.2 0.6 0.2; 0.3 0 0.7];

% chain length
chain_length = 100;

% number of individuals
myN=10000;

% containers
chain = zeros(1,chain_length);
a1=chain*0.;
na1=a1;
a2=chain*0.;
na2=a1;
a3=chain*0.;
na3=a1;

mymean=2;

for j=1:myN
     % constructing markov chaing
     chain = zeros(1,chain_length);
     chain(1)=(randn()<0.5)+(randn()<0.5)+1;
     for i=2:chain_length         
          this_step_distribution = transition_probabilities(chain(i-1),:);
          cumulative_distribution = cumsum(this_step_distribution);
          r = rand();
          chain(i) = find(cumulative_distribution>r,1);
     end

     % the mean for generating fixed effects
     mymean=mymean+3*randn();

     % generating fixed effects  
     fe=mymean+randn();
%      disp(fe)

     % counting individuals in each state
     na1=na1+(chain==1);
     na2=na2+(chain==2);
     na3=na3+(chain==3);

     % summing the individuals
     a1(chain==1)=a1(chain==1)+fe;
     a2(chain==2)=a2(chain==2)+fe;
     a3(chain==3)=a3(chain==3)+fe;
end

% constructing means in each group
a1=a1./na1;
a2=a2./na2;
a3=a3./na3;

% calculating weights
w1=na1/myN;
w2=na2/myN;
w3=na3/myN;

% weighting
wa1=w1.*a1;
wa2=w2.*a2;
wa3=w3.*a3;

% check for weights to sum
w1+w2+w3;

% check that the aggregate is fixed across time
wa1+wa2;
wa1+wa2+wa3;

% checking the relationship

% estimates of mc probabilities for 2 states
[wa1(1:end-1); wa2(1:end-1)]'\wa1(2:end)'
% estimates of mc probabilities for 3 states
[wa1(1:end-1); wa2(1:end-1); wa3(1:end-1)]'\ wa1(2:end)'



% estimates of mc probabilities for 3 states
[w1(1:end-1); w2(1:end-1); w3(1:end-1)]'\ w1(2:end)'

[a1(1:end-1); a2(1:end-1); a3(1:end-1)]'\ a1(2:end)'
[a1(1:end-1); a2(1:end-1); a3(1:end-1)]'\ a2(2:end)'
[a1(1:end-1); a2(1:end-1); a3(1:end-1)]'\ a3(2:end)'