Yes, I was going to suggest you work on the double expansion for what you need. That can be automated to very high order. In general, the "exact" terms [$]a_i(X,\sigma_0)[$] are very difficult to find analytically, while in contrast the double expansion terms can be found analytically with a good computer algebra program like Mathematica.
I will add that, in the abstract, the exact expansion certainly 'works' in the sense that the power series exists. It is a general consequence of the so-called heat kernel expansion for n-dimensional diffusions. But, as I said, actually getting an analytic expression, beyond say [$]O(T)[$], can be very difficult. Just look at what Paulot
had to do to get the SABR case done exactly through [$]O(T^2)[$] and you'll see the difficulties. Likely some numerics with geodesics would get you the terms [$]a_i(X,\sigma_0)[$], numerically. But, if you're going to go that route, you might as well just solve the full evolution PDE numerically.