Curvature and Parallel Transport Here is an updated formulation of the question, which is more precise and I think completely correct:
Suppose $M$ is a Riemannian manifold. Pick a point $p$ in $M$ and let $U$ be a neighborhood of the origin in $T_p M$ on which $exp_p$ restricts to a diffeomorphism. Let $X$ and $Y$ be tangent vectors in $T_p M$, and let $V$ be the intersection of $U$ and the plane spanned by $X$ and $Y$. Let $c(t)$ be a piecewise smooth simple closed curve in $V$. I claim that for any vector $Z$ in $T_p M$,
$R(X,Y)Z=(P_c(Z)−Z)Area(c)+o(Area(c))$
where $R$ is the Riemannian curvature tensor, $P_c$ is parallel transport around the image of $c$ under $exp_p$, and $Area(c)$ is the area enclosed by the image of $c$ under $exp_p$.
Can anyone refer me to a proof of this statement or something similar? I am fairly sure the argument has something to do with integrating the curvature 2-form over the embedded surface obtained by restricting $exp_p$ to the region enclosed by $c$, but I am having trouble with the estimates. Unfortunately I can't find anything in Kobayashi and Nimazu.
Thanks in advance!
Paul
 A: It appears to me that one reason why nobody has proved the formula yet is that the formula is still wrong. First, the formula has to depend on $X$ and $Y$. If you rescale $X$ and $Y$, the left side of the formula scales but the right side stays constant. That can't be. Second, the two sides of the equation do not scale the same under a constant scaling of the metric.
Warning: I wrote this up very quickly and did not check for typos and errors. It's possible that my final formula is still not right, but I am confident that my argument can be used to obtain a correct formula. I also did not provide every last detail, so, if you're unfamiliar with an argument like this, you need to do a lot of work making sure that everything really works. The key trick is pulling everything back to the unit square, where elementary calculus can be used. I'm sure this trick can be replaced by Stokes' theorem on the manifold itself, but that's too sophisticated for my taste.
Holonomy calculation
ADDED:
The correct formula, if you assume $|X\wedge Y| = 1$, is
$P_\gamma Z - Z = Area(c) R(X,Y)Z$
This scales properly when you rescale the metric by a constant factor. Notice that the left side is invariant under rescaling of the metric.
I recommend looking at papers written by Hermann Karcher, especially the one with Jost on almost linear functions, the one with Heintze on a generalized comparison theorem, and the one on the Riemannian center of mass. I haven't looked at this or anything else in a long time, but I have the impression that I learned a lot about how to work with Jacobi fields and Riemann curvature from these papers.
Finally, don't worry about citing anything I've said or wrote. Just write up your own proof of whatever you need. If it happens to look very similar to what I wrote, that's OK. I consider all of this "standard stuff" that any good Riemannian geometer knows, even if they would say it differently from me.
EVEN MORE: There are similar calculations in my paper with Penny Smith: P. D. Smith and Deane Yang
Removing Point Singularities of Riemannian Manifolds, TAMS (333) 203-219, especially in section 7 titled "Radially parallel vector fields". In section 5, we attribute our approach to H. Karcher and cite specific references.
A: These matters are discussed in great detail in lecture 19 of the book "Differential geometry" by Postnikov. The original Russian edition is reviewed in [MR0985587 (90h:53002)]; there is also French edition, but I am not aware of an English one.
A: Your formula for $R(X,Y)Z$ seems to be identical to the formula on page 256 of Peter Petersen, "Riemannian Geometry", second edition, Springer 2006. Petersen gives a sketch of a proof, and calls the formula "Cartan's characterization of the curvature".
A: I believe, this is called Ambrose-Singer theorem. For the proof - you may introduce some coordinates (s,t) in the exp-image of plane spanned by X and Y, and define V(s,t) to be parallel along, say s-coordinate lines. Then compute how the derivative of V in t-direction changes along s-coordinate lines: it is D_s D_t V = R(X,Y)V since D_t D_s V \equiv 9, and s, t coordinate vectors commute - then integral of D_t V (s,t) = \int R(X,Y)V - D_t V(0,t) which is your parallel transport ... it may be something about this in doCarmo Riemannian geom, or Milnor Morse theory ...
A: I would put this in a comment if I had enough reputation points, but something's wacky about your formula. From the definitions you gave, the RHS has no dependence on X and Y explicitly, but the LHS is tensorial in X, Y. If you replace X and Y by 2X and -Y, the plane spanned by them will be the same, and so the domain V is the same. So the RHS doesn't change and the LHS becomes -2 of what it is. 
While I can see the problem with the sign-change associated to swapping X and Y be taken care of if you used the signed area, I think you need to clarify some of your definitions for the formula to make sense. 
To be more precise, the second term on the RHS is small o, so we will ignore it. The first term is the area times a holonomy element. By definition of parallel transport, $\|P_cZ - Z\| \leq 2 \|Z\|$ by the triangle inequality, so for a fixed Z the RHS is $\lesssim Area(c)$, which if you can make arbitrarily small. The left hand side for fixed $X,Y$ is, well, fixed. 
