restart:

# Modelling
# Let A, B, C, D, E the symbols that represent each family and let us consider
# each of them has a bow witin wich will be put one of these 5 letters to indicate
# the origine of the gift the owner of the box is boing to receive.
#
#      A1 A2 A3   B1 B2 B3 B4   C1 C2 C3 C4   D1 D2 D3 D4   E1
#
# If one is only interested in computing the probability of, let us say,
# all the members of family B receive gifts from members of faniliy C, one
# can forget the numbers 1, 2, 3, 4 and simply consider these boxes:
#
#      A A A   B B B B   C C C C   D D D D   E
#
# The probability P the whole family B will receive gifts from family C
# means that the content of these 16 boxes are
#
#      x x x   C C C C   x x x x   x x x x   x
#
# where x denotes "anything but C"
#
# There are 16! / (3! * (4!)^3 *1!) to dispose the 16 gofts withib the 16 boxes.
# On the other side there are 4! ways to place the 4 C's gift within the 4 B's boxes
# and 12! waus to place the 12 remaining gifts within the 12 remaining boxes.
#
# Thus the probability that the 4 members of family B receive the 4 C's gift is :

P1 := 4! * 12! / 16!;
evalf(%);

# What is the probability that all members of family B receive gits from family C
# and that all members of family C receive gifts from family B?
#
# Here the picture of the contentd of the boxes is
#
#      x x x   C C C C   B B B B   x x x x   x
#
# where x denotes "anything but B or C".
#
# As before there are 4! ways to put the C's gift within the B's boxes.
# There are also 4! ways to put the B's gift within the C's boxes.
#
# Let us consider the deposit of C's gifts within B's boxes.
# The 8 remaining gifts can be placed in 8! different ways in the boxes
# labeled x. Once done, the 4 B's gift can only be places in C's boxes in
# 4! different ways.

P2 := 4! * 4! * 8! / 16!;
evalf(%);

# Note that these two results can be directly obtained this way
# (but reasoning cannot harm):

P1 := 1 / combinat:-multinomial(16, 4, 12);
P2 := 1 / combinat:-multinomial(16, 4, 4, 8);

# Finally let us go a litle bit further than problem 2 by imposing now
# that if gift Ci is put within box Bj, then gift Bj is put within box Ci.
#
# Let us start from the result of Problem 2
# We got P2 := 4! * 4! * 8! / 16! because they were 4! to place C's gifts
# in B's boxes and 4! to place B's gifts in C's boxes.
# But now the constrain "if box Bi contains gift Cj then box Cj contains
# gift Bi" reduces this 4! * 4! number.
#
# In effect, once an arrangement is done for, let us say, the content of
# each B box, then the arrangement of the gifts in each box is a direct
# consequence of the former.
# For instance
#      if B1 has C4, B2 has C3, B3 has C1, B4 has C2
# which can be written S := [1=4, 2=3, 3=1, 4=2],
# then the content of C is defined by (rhs=lhs)~(S) in Maple's language.
#
# Thus

P3 := 4! * 8! / 16!;
evalf(%);

restart:

with(Statistics):

# I describe the families by indexing their members for clarity purpose,
# but this is completely useless for the problem posed.

F  := [seq(A[i], i=1..3), seq(B[i], i=1..4), seq(C[i], i=1..4), seq(D[i], i=1..4), E[1]];
N  := numelems(F);

# Estimation of the number of occcurrences of EV in a random assignement.
#
# There are about 22% of chanced that EV doesn't occur

R := 10^5:
Coincidences := Vector(R):

for n from 1 to R do
  P  := Matrix(8, 2, combinat:-randperm(F));
  M1 := map2(op, 0, P[.., 1]):
  M2 := map2(op, 0, P[.., 2]):
  is~(M1=~M2);
  Coincidences[n] := numelems(Select((x -> is(x=true)), %))
end do:

sort(Tally(Coincidences), key=(x -> lhs(x)));
Histogram(Coincidences, discrete=true, thickness=11)

# Results
#
# # There are about 22% of chanced that EV doesn't occur

# Estimation of the number of occcurrences of EV in a random assignement.

R := 10^5:
IntraExchanges := Vector(R):

for n from 1 to R do
  P  := Matrix(8, 2, combinat:-randperm(F));
  P0 := map2(op, 0, P):
  L0 := map(sort, convert(P0, listlist));
  IntraExchanges[n] := ListTools:-Occurrences(2, [seq(ListTools:-Occurrences([f, f], L0), f in [B, C, D])]);
end do:

sort(Tally(IntraExchanges), key=(x -> lhs(x)));

# Results
#
# 4373 times out of 10^5 there is  1 family   for which all the members exchange with themselves.
#  109 times out of 10^5 there are 2 families for which all the members exchange with themselves.
#    3 times out of 10^5 there are 3 families for which all the members exchange with themselves.

# Estimation of the number of occcurrences of EV in a random assignement.

ff := [[B, C], [B, D], [C, D]]:
R  := 10^5:
EV := 0:

for n from 1 to R do
  P  := Matrix(8, 2, combinat:-randperm(F));
  P0 := map2(op, 0, P):
  L0 := map(sort, convert(P0, listlist));
  [seq(ListTools:-Occurrences(f, L0), f in ff)];
  if has(%, 4) then EV := EV+1 end if:
end do:

EV

# Results
#
# The frequency of event EV is about 4e-3.
# (Note that the frequency of EV for a given couple of families is 4e-3/3).
#
# Assuming EV happened last Christmas between families, let's say B anc C, what is the
# (estimation of) that EV happens this Christmas between the same families again?
#
# If EV happens this Christmas it can concern families (B, C) or (B, D) or (C,D).
# The probability that EV happens and concerns families B and C is about
# 4e-3 * 1/3.
#
# The probability that EV appears twice in a row between families B and C is
# (4e-3 / 3)^2 ~ 1.8e-6
#
# Quite small a value.
# So, if I understood correctly the problem, the probability that all the members of
# family B exchanged with all the members of family C twice in a row is a strong
# evidence that something not random did occur.