These are the basic steps you need to run through (is this what you have done? Apologies if this is too basic)...
Make sure that your 1940 shapefile has an attribute containing county
area (e.g. ‘Area1940’)
use the Union tool (ArcToolbox > Analysis Tools > Overlay > Union)
with your 1880 and 1940 shapefiles as inputs.
Specify a new output shapefile.
In the 'JoinAttributes' option select 'ALL'.
Leave the XY tolerance option blank (see later).
Check 'Gaps allowed'
This should produce an output shapefile in which each polygon has attributes from both the 1880 and 1940 shapefiles.
You will need to calculate the areas of your new polygons, so open the table and add a new field (of 'double' type), and use the field > calculate geometry to calculate the new areas (e.g. ‘ AreaUni’).
You can then add another field in which you can calculate percentages, which should be: (AreaUni/Area1940)*100
It might help to decide exactly which fields you will need in your output shapefile before you do the union. How many fields do your two input shapefiles have? For this analysis you really only need county names or Ids from each time and area for 1940. It might make things clearer if you make copies of your shapefiles and delete any unnecessary fields from them before doing the union.
But there are several issues which may complicate things. The biggest is probably how the two boundary sets were digitised. If boundaries that should be the same are slightly different you will end up with hundreds of small slither polygons along the edges of your union output. These can be dealt with during the union process by specifying an XY tolerance. Specifying a tolerance will snap together verticies that are within the tolerance, thereby eliminating slither polygons. I’m not sure that there is a guide for what value to use – you might have to do a bit of trial and error, and have a look at the differences between like boundaries.