Understanding $O_p$ One thing I feel like I have never mastered is the concept of $O_p$ convergence and how to use it. I understand the basic idea and what bounded in probability means, but I always have a hard time understanding how to apply it myself. I have an example which I hope someone can explain to me.
Actually, I have two related examples.  If we let $\epsilon_t$ be $N(0, I)$, consider:
$$
A=\frac{1}{T}\sum_{t=1}^T\epsilon_t\sum_{i=1}^{t-1}\epsilon_i'\\
B=\frac{1}{T}\sum_{t=1}^T\sum_{i=1}^{t-1}\epsilon_i\sum_{j=1}^{t-1}\epsilon_j'
$$
I know that $A\in O_p(1)$ and that $B\in O_p(T)$, but I want to fully understand why. Can anyone explain why they are of these orders? Preferably in an intuitive rather than a  rigorous way. 

As for my own thoughts, for $A$, the way I'm thinking is that $\epsilon_t\sum_{i=1}^{t-1}\epsilon_i'$ is $O_p(1)$, so essentially what we have is $T$ $O_p(1)$ terms. Since we also divide by $T$, what remains is $O_p(1)$. 
For $B$, I'm not sure how I can see that this is true. One idea is that we have two terms (the two $i$ and $j$ sums) that increase linearly with $T$, so thus we get $O_p(T^2)$, which is divided by $T$ to give $O_p(T)$. But is it true in general that $O_p(T)O_p(T)=O_p(T^2)$?

Edit: For my last question, if it is true in general that $O_p(T)O_p(T)=O_p(T^2)$, I think it is. If $X_T\in O_p(T)$, then that means $X_T=Y_T\times T$, where $Y_T\in O_p(1)$. So then
$$
O_p(T)O_p(T)=TO_p(1)TO_p(1)=T^2O_p(1)=O_p(T^2).
$$
 A: To get the technicalities over with
$X_n = \mathcal{O}_p(a_n)=\text{Pr}\left(X_n/a_n > M\right) < \delta$ for every $\delta > 0$.  I'll also assume that the sequence is i.i.d.
Now the way I approach both of these problems is graphically.  The first problem is really an increasing triangular array.
$\begin{array}{lllll}
  & 1 & 2 & \dots & T\\
1 & \epsilon_1\epsilon_1  &   &       &  \\
2 & \epsilon_2\epsilon_1  & \epsilon_2 \epsilon_2   &       & \\
 \vdots &  &   &  \ddots     & \\
T & \epsilon_T\epsilon_1  & \epsilon_T\epsilon_2 & \dots & \epsilon_T\epsilon_T
\end{array}$ 
The expectation of this array is simply
$\begin{array}{lllll}
  & 1 & 2 & \dots & T\\
1 & \sigma_{11}  &   &       &  \\
2 & 0  & \sigma_{22}   &       & \\
 \vdots &  &   &  \ddots     & \\
T & 0  & 0 & \dots & \sigma_{TT}
\end{array}$ 
if each of the $\sigma_{ii}$ are bounded then the sum of the expectations divided by $T$ are bounded by the max $\frac{T \sigma_{ii}}{T}$ times some finite constant that can be adjusted to meet any desired $\delta > 0$. 
The second problem is quite similar, except for each $i$ you are adding triangular arrays of size $i$.  This means that you have $T(T+1)/2$ non-zero expectations.  For $\mathcal{O}_p$ like other big-O notation only the highest order polynomial matters, so the sum of $B$ is bounded by $T$.
If this is a bit confusing, try to sit down with pen and paper or excel and work out a simple example with T = 3.  This should help a lot with your understanding.
For your Edit, recall that for the last statement to be true 
$\int_{-\infty}^{M_1T} dX_n \int_{-\infty}^{M_2T}  dY_n = \int_{-\infty}^{M_3T^2} d(X_n Y_n)$ for some integers $M_1$, $M_2$, and $M_3$.
This statement being true depends on how you define $Y_n$ and $X_n$, if they are independent independent sequences e.g. any element of $X_n$ is independent of every other element of $X_n$ and $Y_n$, then the probability distribution can be broken down into the product of the marginal distributions and you are good-to-go.  
If you drop normality though things also start falling apart a $T$ distribution with 2 degrees of freedom would be bounded in probability, but it's product would not be.
A: If $A_T \in \mathcal{O}_p(1)$, that means $\{A_T/1\}$ is uniformly tight, or bounded in probability. Notice that I added a $T$ subscript. That's because this is a property for a family of random variables. The family in our cases is $A_2, A_3, \ldots$. This property is saying that all of these random variables are "pretty much" bounded above and below by the same number. By "pretty much" I mean that they are bounded probabilistically, meaning, the probability of breaking out of these bounds is negligible. We can increase that bound in order to make the probability smaller than whatever we want.
Notice I'm changing the bounds of the outer sum, here. $t$ needs to start at $2$, or else the inner sum doesn't make any sense. Also, I removed the "prime" symbol from the inner $\epsilon$s, because it is unnecessary.
The good news is that you don't need triangular arrays to show the first one. Let $\epsilon > 0$, then we can show there exists an $M > 0$ such that 
$$
P(|A_T| \ge M) < \epsilon
$$
for all $T$. Just pick $M > \epsilon^{-2}$ and notice that
\begin{align*}
P(|A_T| \ge M) &= P\left(\left|\frac{1}{T}\sum_{t=2}^T\epsilon_t\sum_{i=1}^{t-1}\epsilon_i\right| \ge M \right) \\
&\le M^{-2}T^{-2}\mathbb{V}\left[\left|\sum_{t=2}^T\epsilon_t\sum_{i=1}^{t-1}\epsilon_i\right| \right] \tag{Chebyshev's} \\
&\le M^{-2}T^{-2}\mathbb{V}\left[\sum_{t=2}^T\epsilon_t\sum_{i=1}^{t-1}\epsilon_i \right] \tag{*} \\
&= M^{-2}T^{-2}\left\{\sum_{t=2}^T \mathbb{V}\left[\epsilon_t\sum_{i=1}^{t-1}\epsilon_i\right] + 2\sum_{t=2}^{T-1}\sum_{s=t+1}^T \mathbb{Cov}\left[\epsilon_t\sum_{i=1}^{t-1}\epsilon_i ,\epsilon_s\sum_{j=1}^{s-1}\epsilon_j \right] \right\} \tag{bilinearity} \\
&= M^{-2}T^{-2}\left\{\sum_{t=2}^T \mathbb{E}\left[\left(\sum_{i=1}^{t-1}\epsilon_i\right)^2\right] + 2\sum_{t=2}^{T-1}\sum_{s=t+1}^T \mathbb{Cov}\left[\epsilon_t\sum_{i=1}^{t-1}\epsilon_i ,\epsilon_s\sum_{j=1}^{s-1}\epsilon_j \right] \right\} \tag{zero mean}\\
&= M^{-2}T^{-2}\left\{\sum_{t=2}^T (t-1) + 2\sum_{t=2}^{T-1}\sum_{s=t+1}^T \mathbb{Cov}\left[\epsilon_t\sum_{i=1}^{t-1}\epsilon_i ,\epsilon_s\sum_{j=1}^{s-1}\epsilon_j \right] \right\} \tag{normality} \\
&= M^{-2}T^{-2}\left\{(T^2+T)/2 -T + 2\sum_{t=2}^{T-1}\sum_{s=t+1}^T \mathbb{Cov}\left[\epsilon_t\sum_{i=1}^{t-1}\epsilon_i ,\epsilon_s\sum_{j=1}^{s-1}\epsilon_j \right] \right\} \tag{algebra} \\
&= M^{-2}T^{-2}\left\{(T^2+T)/2 -T + 2\sum_{t=2}^{T-1}\sum_{s=t+1}^T \mathbb{Cov}\left[\sum_{i=1}^{t-1}\epsilon_t\epsilon_i , \left\{\sum_{j=1}^{t-1}\epsilon_s\epsilon_j + \epsilon_s\epsilon_t\right\} \right] \right\}\\
&= M^{-2}T^{-2}\left\{(T^2+T)/2 -T  \right\} \\
&\le M^{-2}
\end{align*}
For the line with (*), I'm using the fact that $\mathbb{V}[|X|] \le \mathbb{V}[X] + (\mathbb{E}[X])^2$. For the line above that, you can use bilinearity, recognize where the means are zero, and then independence. 
Regarding 
$$
B=\frac{1}{T}\sum_{t=2}^T\sum_{i=1}^{t-1}\epsilon_i\sum_{j=1}^{t-1}\epsilon_j' = \frac{1}{T}\sum_{t=2}^T (t-1) z_t z_t'
$$
you can use Chebyshev's on that again.
