How can I adjust values in a timeseries to account for effects from other variable(s)? i.e. using a GAM I have a dataset of monthly air pollution levels (i.e. PM2.5) as well as a corresponding dataset with monthly rainfall amount (i.e mm) spanning several years. These datasets have moderately strong, negative relationship (i.e. R-squared = 0.6). I show an example of the data below.
I have read in several articles where people describe adjusting ozone levels based on meteorology effects using Generalized Additive Models (GAM) (such as this article) 
However, I don't really understand how this is done as the articles are not very descriptive. For example, which functions and parameters would be selected in R or python GAM implementations and why? Usually the examples I find only provide examples of smoothing.
Any thoughts or suggestions or tutorial links would be greatly appreciated. I am also open to other suggestions which use different methods aside from GAMs.
 A: Since you asked me to work my special magic on your data using automatic software ... I introduced your 60 monthly values for two time series to AUTOBOX , my preferred tool of choice ( since I helped develop it). No hint whatsoever was provided to the system regarding a theoretical expectation of seasonal factors in BC .
Adjusting one series for 1 or more supportive series is basically the predicted value from a transfer function (regression + ). Predicted values need to be sometimes adjusted for exogenous (unspecified outside) effects whose form needs to be established after adJusting for memory (ARIMA). If you wish to post your data I will try to  help further.
EDITED AFTER RECEIPT OF DATA :
The idea is simple enough ... form a useful model between BC and Rainfall.  Following https://web.archive.org/web/20160216193539/https://onlinecourses.science.psu.edu/stat510/node/75/ while being concerned with possible anomalies and changes in model error variance http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html we use the following flow diagram as a paradigm http://autobox.com/cms/images/dllupdate/TFFLOW.png for the steps in the analysis.
The transfer function equation suggests non-constant error variance and a few anomalies pulses and 1 level shift and three very definite deterministic/fixed seasonal effects.   i.e. 3 seasonal pulses (march/june/november) . The original plot of Y (BC) is here  and X (RAINFALL) is here  . Using Tsay's method if testing for constant error process variability ( always a good idea ! ) we obtained  time period 48 as a break point given the effect of X and the three seasonal dummies.
A plot of the residuals form the model ( another good idea ! ) suggest sufficiency  as does the ACF of the errors  . The ACTUAL and the CLEANSED DATA is shown  here   . At the end of the day the FITTED VALUES from the model shown here with the ACTUAL  and partially here in tabular form  reflect BC adjusted for RAINFALL and the identified repetitive (seasonal pulse; level shift ) and one-time only deterministic components (pulses).
In effect we have identified the need to adjust BC for Rain while also adjusting for seasonality and anomalous values & increased variability over time starting at period 48. The automatic identification  of the three seasonal dummies supports your hypothesis thus all is good ... data analysis can be used to support prior degrees of belief ! These seasonal factors represent ADDITIONAL IMPACT above and beyond the normal response of BC to RAINFALL.
