The formula for nn.subsampling looks like this:
output[i][j][k] = bias[k] + weight[k] * sum_{s=1}^kW sum_{t=1}^kH input[dW*(i-1)+s)][dH*(j-1)+t][k]

We can rearrange the formula in order to use some of the existing modules:
O1[i][j][k] = 1/(kH * kW) * sum_{s=1}^kW sum_{t=1}^kH input[dW*(i-1)+s][dH*(j-1)+t][k] # This is AvgPool2d
O2[i][j][k] = O1[i][j][k] * (kH * kW) # Multiply the matrix by a constant
O3[i][j][k] = O2[i][j][k] * weight[k] + bias[k] # Use a learnable vector `weight` and `bias` to calculate O3