smat = {};
For[j = 1, j < 3, j++,

  ClearAll[zvec, UVect, uvec];
  zvec[-4] = ConstantArray[0, 12];

  sigmaumat3 = IdentityMatrix[20];
  UVect[-4] =
   RandomVariate[
    MultinormalDistribution[ConstantArray[0, 20], sigmaumat3]];
  uvec[t_] :=
   uvec[t] =
    RandomVariate[
     MultinormalDistribution[ConstantArray[0, 4], IdentityMatrix[4]]];
  UVect[t_] := UVect[t] = Drop[Join[uvec[t], UVect[t - 1]], -4];
  zvec[t_] := zvec[t] = azmat.zvec[t - 1] + bumat.uvec[t];
  Zmat[t_] := Zmat[t] = Join[UVect[t], zvec[t - 3]];

  smat = Join[smat, {Zmat[j]}];
  ];
smat // MatrixForm