Prove the AGM identity using only Hypergeometric series Arithmetic Geometric Mean can be represented by a Hypergeometric function:
$$\text{agm}(1,p)=\frac{1}{{_2F_1} \left(\frac{1}{2},\frac{1}{2};1;1-p^2 \right)}$$
$$0<p \leq 1$$
One of the main properties of the AGM is the following identity:
$$\text{agm}(1,p)=\frac{1+p}{2}\text{agm} \left(1,\frac{2\sqrt{p}}{1+p} \right)$$
This allows the infinite product representation of the AGM.

I wanted to know if it's possible to prove this identity by directly using the Hypergeometric series.

For the Hypergeometric function the identity will take the following form:
$${_2F_1} \left(\frac{1}{2},\frac{1}{2};1;1-p^2 \right)=\frac{2}{1+p} {_2F_1} \left(\frac{1}{2},\frac{1}{2};1;\frac{(1-p)^2}{(1+p)^2} \right)$$
In the series form it will be:
$$ \sum_{k=0}^\infty \frac{1}{k!^2} \left(\frac{1}{2}\right)^2_k (1-p^2)^k=\frac{2}{1+p} \sum_{k=0}^\infty \frac{1}{k!^2} \left(\frac{1}{2}\right)^2_k \frac{(1-p)^{2k}}{(1+p)^{2k}}$$
I think the following substitution will simplify things:
$$1-p=2x$$

$$ \sum_{k=0}^\infty \frac{1}{k!^2} \left(\frac{1}{2}\right)^2_k 2^{2k}x^k(1-x)^k=\frac{1}{1-x} \sum_{k=0}^\infty \frac{1}{k!^2} {\left(\frac{1}{2}\right)^2_k} \frac{x^{2k}}{(1-x)^{2k}}$$

I haven't been able to prove this identity from the series.
Comparing terms in this form is useless, since the partial sums of the series are not equal (the second series converges much faster).
The only idea I have is to use the uniqueness of the power series, which requires expanding everything, so there are only powers of $x$ left.
We have:
$$(1-x)^k=\sum_{l=0}^k (-1)^l \left(\begin{matrix} k \\ l \end{matrix} \right) x^l$$
$$\frac{1}{(1-x)^{2k+1}}=(2k)! \sum_{n=0}^\infty (-1)^{n+2k}~ (n+2k)_{2k} ~x^n$$
Here $(n+2k)_{2k}$ actually means falling factorial, not rising factorial, like above. $(n+2k)_{2k}=(n+2k)(n+2k-1)(n+2k-2) \cdots$. I don't know what other notation to use in this case.
Now I'm stuck. I don't know how to get the single power series for $x$ on each side so we can compare them.

I think it is faster to go through an elliptic integral. Given $a,b\in\mathbb{R}^+$ we may define
$$ E(a,b)=\int_{0}^{+\infty}\frac{dx}{\sqrt{(a^2+x^2)(b^2+x^2)}} \tag{1} $$
and notice that by Lagrange's identity $(a^2+x^2)(b^2+x^2)=(ax+bx)^2+(ab-x^2)^2$ and a suitable change of variable we have
$$ E(a,b)=E\left(\frac{a+b}{2},\sqrt{ab}\right) \tag{2}$$
hence:
$$ E(a,b)=E(AGM(a,b),AGM(a,b))=\frac{1}{AGM(a,b)}\int_{0}^{+\infty}\frac{dx}{x^2+1}\tag{3} $$
and 
$$ AGM(a,b) = \frac{\pi}{2 E(a,b)}\tag{4} $$
so our claim is equivalent to proving that 
$$\frac{2}{\pi}\int_{0}^{+\infty}\frac{dx}{\sqrt{(1+x^2)(p^2+x^2)}} = \phantom{}_2 F_1\left(\frac{1}{2},\frac{1}{2};1;1-p^2\right) \tag{5} $$
but by setting $u=\frac{1}{1+x^2}$, that simply follows from Euler's integral representation for $\phantom{}_2 F_1$.
