Introduction
Team GO Rocket’s (TGR) takeover of Pokemon GO provides a new domain of PvE experience using PvP mechanics (contrast to Raids, which use strictly PvE mechanics). Trainers fight in 3v3 combat against TGR Grunts, who use Shadow Pokemon with greatly inflated stats to offset the natural advantage that a human player has.
With CPs in the thousands, the question since day 1 has been: what are the stats of TGR’s Shadow Pokemon? Raid Boss stats were simple to crack, but Shadow Pokemon stats have proven difficult. Previous attempts have yielded progressively better approximations, but no formula that predicts Shadow Pokemon CP exactly. Due to imprecision in floating point numbers and high CP products of Shadow Pokemon, this was attributed to rounding errors compounded in several different places.
This article presents a formula that predicts exactly a TGR Shadow Pokemon’s stats (to an extreme degree of precision) and its derivation.
Assumptions
As there are infinite individual combinations of stats that can produce a CP value, I required a few assumptions to help guide refinement of the formula:

One formula applies to all Shadow Pokemon; that is, individual Shadow Pokemon won’t have random IVs or stat changes.

Applying Occam’s razor, any derived coefficients should be fairly simple, elegant looking numbers and not like, irrational.

When I use multipliers, they are in relation to the species’s stats at L40 (CpM = 0.7903).

Shadow Pokemon, like Raid Bosses, use 15 IVs in all stats.
Methods: Finding Stats
Narrowing Attack range
All other variables known, one can rewrite the damage formula to solve for a range of Attack values for a Shadow Pokemon:
[Min_Atk, Max_Atk) = [Dmg  1, Dmg) * 2 * Def / (1.3 * Pow * STAB * Eff)
As damage sustained is not healed after a TGR Grunt battle, the easiest way to get a damage value is to enter a battle and run away after sustaining 12 Fast Moves. This will yield an exact value to use in the equation. Taking more damage in a single hit will produce a smaller Attack range. Assuming that all Shadow Pokemon use the same Attack multiplier, finding Attack ranges across multiple matchups will narrow down the Attack multiplier further.
Across the small sample size that I had, the Attack multiplier was over the range (3.185, 3.214) compared to a regular L40 15 Atk Pokemon. celandro collected a more extensive data set that settled on an Attack multiplier tightly distributed around 3.201.
Finding the exact HP
In order to figure out the HP and Defense multipliers, we have to know one of those stats first. Otherwise, 2x HP looks almost exactly like 2x Defense, for example. There are a few creative ways to do this.
Yawn to death / yellow HP
This is the most basic way. Yawn and Splash always do 1 HP damage, so KO’ing a Shadow Pokemon with x Yawns means they have x HP. Simple, right? But Yawn and Splash are so bad that either your 3 Pokemon will wipe out before KOing the enemy or the battle will just time out. The only exception was Duskull, who has low HP and was double resisted by Slaking and Snorlax. This was an critical data point, because it lets me use other methods to find HP for all other Shadow Pokemon species.
Alternatively, reducing a Pokemon’s HP bar to yellow means that you’ve dealt >50% damage. This narrows down their HP to 2 possible values.
Division with remainders
For the nonDuskull Pokemon who can’t be KO’d outright by Yawn, you can whittle them down with Fast Moves from another Pokemon until they’re about to be KO’d, then do the remaining damage with Yawn. You won’t know how much damage the other Pokemon is doing, but that doesn’t matter.
If you know that Ralts is KO’d by 8 Tyranitar Bites + 6 Yawns or by 3 Kyogre Waterfalls + 1 Yawn, then it must have 70, 94, 118, etc. HP. Since we know from the above example that Duskull has 61 HP, and assuming that Ralts has a similar ratio to its (base HP + 15), then Ralts must have 70 HP. Even without knowing exactly how much damage each move does, if you take into account stats and move power, then 70 HP is the only value that makes sense.
Doing 50% HP exactly
In a PvP battle, when a Pokemon is at 50% HP exactly, its HP bar stays green. So if Shuppet’s HP bar is green after 4 Tyranitar Bites, but it faints to 8 Tyranitar Bites, then Shuppet’s HP must be exactly divisible by 8. Assuming that Shuppet has the similar ratio to its (base HP + 15), then Shuppet must have 88 HP.
Species  Base HP  Actual HP  Min. Multiplier  Max Multiplier 

Duskull  85  61  0.772  0.785 
Abra  93  67  0.785  0.797 
Ralts  99  70  0.777  0.788 
Zubat  120  84  0.787  0.797 
Shuppet  127  88  0.784  0.793 
Larvitar  137  94  0.783  0.791 
Finding exact HP values was quite time consuming, so I only collected a small sample. If you look closely, you’ll notice that the minimum HP multipliers for some species is larger than the maximum HP multiplier for Duskull. So something’s going on, maybe one of those rounding errors.
Narrowing Defense range
All other variables known, one can rewrite the damage formula to solve for a range of Defense values for a Shadow Pokemon:
(Min_Def, Max_Def] = 1.3 * Pow * STAB * Eff * Atk / (2 * (Dmg, Dmg  1])
Once you know a Shadow Pokemon’s HP, you can figure out how much damage any Fast Move does by using the division with remainders method above. For example, if Larvitar is KO’d by 4 Feraligatr Waterfalls + 1 Yawn, then Feraligatr’s Waterfall must do 34 HP damage. As with Attack, dealing more damage in a single hit will produce a smaller Defense range, and finding Defense ranges across multiple matchups will narrow down the Defense multiplier.
Across the small sample size that I had, the Defense multiplier was over the range (1.593, 1.602) compared to a regular L40 15 Def Pokemon.
Methods: Finding a Fit
So now that we have fairly narrow multiplier ranges for HP, Attack, and Defense, it’s time to plug those into the CP formula and see if it works out. We’ll start with an HP mult = 0.7853, Atk mult = 3.201, Def mult = 1.602 (I didn’t really start out with this, but this is to avoid wasting your time).
So we have exact matches for a few species, but everything else is slightly off, some by a little, some by a lot. Also, the predicted HP doesn’t exactly match the actual HP. I was ready to throw in the towel after playing with the multipliers a bit. Is this a rounding error? See if you can notice some details about the data before reading further.
Observation 1: odd / even base HP discrepancy
With the present multipliers, the formula tends to underestimate CP. If you look really really closely, it underestimates species with even base HP more than does species with odd base HP (compare Rattata and Zubat to Abra and Dratini). This means that somehow, even base HP species have slightly more HP than expected, or odd base HP species have slightly less HP than expected.
It turns out that if you round down all odd base HPs to the next even base HP and increase the HP multiplier accordingly, much of this discrepancy goes away. This is using HP mult = 0.7925, with the table sorted by base HP.
Observation 2: higher base HP discrepancy, and
Observation 3: disproportionate stats discrepancy
This looks better, but we aren’t even close to a perfect fit, and there are obvious problems at the bottom of the table where the high base HP species are. The current formula is underestimating their HP. Also, among species with the same base HP, the formula is underestimating those who have relatively higher stats in Atk and Def (compare Hitmonlee and Hitmonchan to Larvitar and Vibrava).
Simply upping the HP multiplier doesn’t fix the problem; then the formula would overestimate HP for species with low base HP and / or relatively less Atk and Def. Because of the nature of this discrepancy, here we have to challenge our assumption that all IVs  especially HP  are 15 by default. We haven’t reconciled the discrepancy between predicted and actual HP either, by the way. Since we have nowhere to go but down, let’s try an HP IV of 14 and bump HP mult to 0.8:
Holy crap! Almost an exact match! And now our HP values correlate, too. Just a little adjustment of the Atk or Def multipliers and we should be in business:
Conclusion
A TGR Grunt’s Shadow Pokemon’s stats can be predicted, to an extreme degree of precision, with the following formulas:
HP = floor(0.7903 * 0.8 * (2 * floor(0.5 * base_HP) + 14))
Atk = 0.7903 * 3.201 * (base_Atk + 15)
Def = 0.7903 * 1.60174 * (base_Def + 15)
CP = floor(0.1 * (Atk * Def^0.5 * (0.7903 * 0.8 * (2 * floor(0.5 * base_HP) + 14))^0.5))
Hold on, we’re not done yet. The HP, Atk, and Def multipliers are nearly in an exact 1:4:2 ratio. Abiding by the principle of Occam’s razor, can we pretty up the coefficients?
Extreme finetuning
(thanks to celandro for his help with this exercise)
It turns out that HP, Atk, and Def multipliers of 0.80034175, 3.201367, and 1.6006835 fit the data perfectly while also being in a 1:4:2 ratio. That’s a lot of decimal places, though. Let’s revisit the assumption that multipliers are in relation to a CpM = 0.7903. CpM, after all, is just a global multiplier to all stats, and Raid Bosses don’t use the same CpMs as trainers.
If we discard CpM = 0.7903, we can get a prettier perfect fit to the data with the following formulas:
HP = floor(0.63251 * (2 * floor(0.5 * base_HP) + 14))
Atk = 0.63251 * 4 * (base_Atk + 15)
Def = 0.63251 * 2 * (base_Def + 15)
CP = floor(0.1 * (Atk * Def^0.5 * (0.63251 * (2 * floor(0.5 * base_HP) + 14))^0.5))
In summary, for L40 trainers, TGR Grunts’ Shadow Pokemon have a global stat multiplier of 0.63251, and then their Attack and Defense and multiplied by 4 and 2, respectively. Their Pokemon hit ~3.2x as hard as a regular L40 Pokemon and have ~28% more bulk.
It is possible that there are still decimal places different from what I’ve proposed. However, any such difference is indistinguishable with the data on hand. These values should be suitable for making very accurate simulations that help us dissect TGR matchups.
A side note
In the process of testing, I tried to calculate Shadow Pokemon’s Attack stats with their Charged Moves, and noticed that Grunts don’t charge their moves fully. Their Charged Move damage is multiplied by around 0.5625  it's as if they hit about half of the circles in the Charged Move minigame!
Special thanks
celandro, for providing a repository of data to greatly narrow down the Attack multiplier
Slaking and Snorlax, for Yawn. I knew maxing a Slaking would come in useful some day.