Influence of Seasonality on Unit root tests I gather from this question and answer, a seasonal series is non-stationary as its mean depends on which month it is.
Suppose I have a series which possibly has a unit root (stochastic trend) but also exhibit strong seasonality (see the chart below).

Now I think I can remove seasonality first and then use a suitable unit root test. My doubt is, would it be theoretically sound if I directly perform unit root test on the original series.
My feeling is that a strong seasonality, such that a large part of variance in series is coming from seasonal fluctuation, can hide the underlying stochastic trend, causing the unit root test to misidentify the series as having no unit root.
Does this make sense? Further, what if we try to handle the seasonality within the unit root test. For example, in ADF test, we include seasonal dummies (just like we include a trend variable to handle deterministic trend)?. But even if this is correct, we wouldn't have the right critical values for it, right?
 A: A few days after I posted this question, I attempted simulations to test my intuition that strong seasonality can hide the underlying stochastic trend in the usual ADF test. Consider the following process:
\begin{align}
y_t &= S_t+y_{t-1} + \varepsilon_t \\
\text{where, } \,\, &S_t \sim N(\mu_i,\sigma_s^2) \,\,\,\text{for }t=4n+i;\tag{$i=1,2,3,4;n \in \mathbb N$}\\
&\sum\limits_i \mu_i = 0; \\
&\varepsilon_t \sim IN(0,\sigma_{\varepsilon}^2)
\end{align}
We can think of this process as a quarterly time series exhibiting stable seasonality. Intuitively, the strongness of seasonality can be measured by $(\text{range}(\mu_i))/\sigma_{\varepsilon}^2$ (further logic in the end). Higher this ratio, more likely it would be that seasonality will hide unit root in the usual test.
Chart below shows the results of the sumulation. 10,000 series of length 1000 are generated from the above model and a pure random walk model. Values for seasonal term are taken as $\mu = (1.7,0.8,-1.0,-1.5)$. DF-statistic is calculated for each series using urca::ur.df(., lag = 0) and density estimates are plotted.

From above chart, it is interesting to see that when $\sigma_{\varepsilon}$ is small the usual unit root test can be very wrong.
Theoretical justification:
In this paper, the authors have derived the analogous DF distribution in presence of an additive outlier. The distribution is similar to the usual DF distribution with an additional term as function of $\theta/\sigma_{\varepsilon}^2$; where $\theta$ is the coefficient of the outlier dummy. The above process is somewhat a special case with four additive type outliers and perhaps that's why we are getting the similar results.
A: I am no expert in this topic, but there seem to be dedicated unit root tests developed for seasonal data.
See for example:
Ghysels, Eric, Hahn S. Lee, and Jaesum Noh. "Testing for unit roots in seasonal time series: some theoretical extensions and a Monte Carlo investigation." Journal of Econometrics 62, No. 2 (1994): 415-442.

While we find the procedure proposed by Hylleberg, Engle, Granger,
and Yoo (1990) the most useful among the alternative procedures,
we caution users of many remaining serious obstacles when testing for
unit roots in seasonal time series.

The Hylleberg, Engle, Granger, and Yoo (1990) approach is also included in the uroot package in R, see here for the Vignette.

According to the ADF test, Dickey et al. (1984) obtained the critical values
of $\phi$ under the following data generating process.
$(1 − L^S)\,y_t = \Phi y_{t−S} + \epsilon_t$, $\epsilon_t ∼ iid (0, \sigma_{\epsilon^2})$
where $S$ is peridocity of the data: 4 in the case of quarterly series and 12
for monthly series. Under the null hypothesis $\Phi$ is equal to 1 and the process
contains all the roots in Table (2), page 7. In practice, however, it may be
some, but not all, seasonal unit roots. Hence, it would be convenient to
specify a model in which the regressors allow to test for individual roots.
That is what Hylleberg et al. (1990) achieve in the case of quarterly series.The seasonal operator can be decomposed into the following polynomials:
$(1 − L^4) = (1 − L)(1 + L)(1 − iL)(1 + iL)$, where the roots are $(±1, ±i)$.
Assuming $y_t$ is generated by an AR(p) process, Hylleberg et al. (1990) show
that $y_t$ can be represented as
$\phi (L) \Delta^4 y_t = \pi_1 y_{1,t−1} + \pi_2 y_{2,t−1} + \pi_3 y_{3,t−2} + \pi_4 y_{3,t−1} + \epsilon_t$
,

Hylleberg, Svend, Robert F. Engle, Clive WJ Granger, and Byung Sam Yoo. "Seasonal integration and cointegration." Journal of Econometrics 44, No. 1-2 (1990): 215-238.
A: A series is not automatically non-stationary if it has a seasonal pattern.
For example, Rob Hyndman shows that a cyclic pattern may exist that makes it unpredictable. There is also a good explanation on towarddatascience.com on that topic.
Instead of directly performing the ADF test I would first consider ACF/PACF plots to see if a single series may be stationary or not, that normally clears out a lot of questions. The correlogram is the first plot to be conducted to get a first glance. Normally this is the starting point in checking if you need additional tests. See here (you have a lot of scrolling to do).
However, I believe there is no harm on using the ADF on a non transformed series, because you should already have insights from your ACF. If you don't have any insights, well I once developed a genetic algorithm which fits on a VAR/VECM in R. These approaches normally need a unified level of variables in the model. Thus if there is one variable that is stationary and one that is not, you should log transform and difference all variables, or difference them a second time, you need a unified level. This is also a sign for me in no harm, as we try to measure before and after logging and differencing with the ADF. In addition, believe in Stata it is also one of the first tests conducted on non transformed variables before the KPSS as this comes after the ADF: If your original variable already has a unit root, then it is non stationary. In the end you want to look at the real data.
Further, what if we try to handle the seasonality within the unit root test. For example, in ADF test, we include seasonal dummies (just like we include a trend variable to handle deterministic trend)?. But even if this is correct, we wouldn't have the right critical values for it, right?
It is possible to do seasonal differencing, search for it with the tag: Rob Hyndman, the master of time series ;-) I do not believe you do it with dummies, you do seasonal differencing and then you can normally interpret ADF.
