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".