Probably not the most efficient way to do so, but for descently sized n, this should work well:

ArrayReshape[
 T[ρ, #] & /@ Tuples[Range[3], n],
 ConstantArray[3, n]
 ]

Probably not the most efficient way to do so, but for descently sized n, this should work well:

T[i_, list_List] := Tr[i.KroneckerProduct @@ PauliMatrix[list]]
n = 4;
ρ = RandomReal[{-1, 1}, 2^n];
A = ArrayReshape[
   T[ρ, #] & /@ Tuples[Range[3], n], 
   ConstantArray[3, n]
   ];
A // Dimensions

{3, 3, 3, 3}