Is the reduced density operator a density operator itself? What this question is about:
Given the situation, that it is only known that a system is in a state $|\psi_i\rangle$ with probability $p_i$, one defines the density operator $$\rho=\sum_i p_i|\psi_i\rangle\langle\psi_i|\tag{1},$$ which describes the state of the system.
If a system $\mathcal{AB}$, composed of two subsystems $\mathcal{A}$ and $\mathcal{B}$, is described by a density operator $\rho$, one defines the reduced density matrix $$\rho_{\mathcal{A}}:=\text{Tr}_{\mathcal{B}}\left(\rho\right)\tag{2}$$ as a description of the state of the subsystem $\mathcal{A}$. This is a useful description of the state of the subsystem because

one is able to calculate the expected value of
$O_{\mathcal{A}}\otimes I$, which is
$\text{Tr}(\rho(O_{\mathcal{A}}\otimes I))$, by calculating the
expecte value of $O_{\mathcal{A}}$, which is
$\text{Tr}(\rho_{\mathcal{A}}O_{\mathcal{A}})$ (which can be shown
using (2)).

The question:
Is the reduced density matrix itself a density matrix? E.g. can it be treated as a density operator of a system just like density operators from the definition (1)? In my opinion this doesn't follow from (2) necessarily, because it's not clear, that the reduced density matrix is (mathematically) an object in the shape of (1).
Since no book explicitly answers this question, I am tempted to believe that the reduced density matrix is a density matrix because of (quote) - so because the calculation of an expected value is analog to the case of density matrices, which implies they are objects of the same "shape".
 A: Density matrices are, by definition, positive semi-definite trace-class operators of unit trace, see for example this. In the following we will only discuss the finite-dimensional case, so we have to show that the reduced density matrix is positive semi-definite and of unit trace.$^\dagger$

To start, consider a complex bipartite Hilbert space $H=H_A\otimes H_B$ of finite dimension. Let $\rho$ be a density operator on $H$. In particular it holds that
$$
 \langle \psi|\rho|\psi\rangle \geq 0 \tag{1}$$
for all $|\psi\rangle \in H$ and
$$\mathrm{Tr} \rho =1 \tag{2} \quad .$$
We define the reduced density operator $\rho_A$ of $\rho$ as
$$ \rho_A:=\mathrm{Tr}_B\,\rho\tag{3} := \sum\limits_{k \in K} \left( \mathbb I_A\otimes\langle \psi_k|\right) \,\rho\, \left(\mathbb I_A\otimes |\psi_k\rangle\right) \quad , $$
for some orthonormal basis $\{|\psi_k\rangle\}_{k\in K}$ in $H_B$, cf. here. Note that $\rho_A$ is a linear operator on $H_A$. We compute its trace:
$$\mathrm{Tr}^{(A)} \rho_A = \sum\limits_{j \in J}\sum\limits_{k \in K} \langle \varphi_j|\otimes \langle \psi_k| \, \rho\, |\varphi_j\rangle \otimes |\psi_k\rangle = \mathrm{Tr} \rho \overset{(2)}{=} 1 \quad , \tag{4} $$
where $\{|\varphi_j\rangle\}_{j\in J}$ and $\{|\varphi_j\rangle \otimes |\psi_k\rangle\}_{j\in J,k\in K}$ are orthonormal bases in $H_A$ and $H$, respectively. We are left to show that $\rho_A$ is positive semi-definite. To do so, define $|\phi_k\rangle := |\varphi\rangle \otimes |\psi_k\rangle$ for an arbitrary $|\varphi\rangle \in H_A$. Then
$$0 \overset{(1)}{\leq}\sum\limits_{k\in K}\langle \phi_k|\rho|\phi_k\rangle = \sum\limits_{k \in K} \langle \varphi| \otimes \langle \psi_k|\,\rho\,|\varphi\rangle \otimes |\psi_k\rangle\overset{(3)}{=}\langle \varphi| \rho_A |\varphi\rangle \tag{5} \quad .  $$
In conclusion, since $\rho_A$ is a (linear) positive semi-definite operator of unit trace on $H_A$, it is a density operator.

$^\dagger$ That this generalizes your definition can be seen as follows (again, working with on a finite-dimensional complex Hilbert space):
If an operator $\sigma$ is positive semi-definite, it is hermitian (cf. the last part of this answer) and thus admits a spectral decomposition:
$$\sigma = \sum\limits_i \lambda_i\, |\lambda_i\rangle \langle \lambda_i| \quad , \tag{6}$$
with $\lambda_i \geq 0$. The trace condition now ensures that $ \sum\limits_i \lambda_i = 1$, so it is of the form given by your equation $(1)$ and hence a density matrix. Conversely, for an operator of the form
$$\sigma := \sum\limits_i p_i\, |\psi_i\rangle \langle \psi_i| \tag{7}$$
with $ \sum\limits_i p_i=1$, $p_i \geq 0$ and $|\psi_i\rangle$ normalized (but not necessarily orthogonal), i.e. a density matrix, it is not hard to prove that $\sigma \geq 0$ and $\mathrm{Tr}\, \sigma =1$.
To summarize: Proving that $\rho_A$ is positive semi-definite and of unit trace shows that it admits a representation of the form $(6)$ and thus proves that it is a density operator in terms of your definition.
