Let us rephrase the question as follows.
Is there a Lagrangian density ${\cal L}$ that (1) implements the global $SO(2)$ symmetry manifestly, and (2) whose Euler-Lagrange equations yield Maxwell's equations in vacuum where there are no sources?
The answer is Yes. Below we will show this. Let the speed of light be $c=1$ from now on.
I) Field variables. The model has $2\times 3=6$ gauge potential fields ${\cal A}^a_i(\vec{x},t)$. Here $i=1,2,3$ are three spatial directions, and $a=1,2$ is an internal $SO(2)$ index. The gauge potential transforms
$${\cal A}^a_i~\longrightarrow~ \sum_{b=1}^2M^a{}_b {\cal A}^b_i \tag{1}$$
in the 2-dimensional fundamental representation of $SO(2)$, where
$$\begin{align}M^a{}_b ~=~\begin{bmatrix} \cos(\theta) & \sin(\theta) \cr -\sin(\theta) &\cos(\theta)\end{bmatrix} ~\in~& SO(2), \cr \sum_{b,c=1}^2 (M^{t})_a{}^b g_{bc} M^c{}_d ~=~& g_{ad}, \cr g_{ab}~\equiv~& \delta_{ab}
.\end{align}\tag{2}$$
The magnetic field $\vec{B}$ and electric field $\vec{E}$ are given by the curl of the gauge potential $\vec{\cal A}^a$,
$$\vec{\cal B}^a := \vec{\nabla} \times \vec{\cal A}^a, \qquad a=1,2, \tag{3}$$
where
$$ \vec{B}\equiv\vec{\cal B}^1 \qquad \mathrm{and} \qquad \vec{E}\equiv \vec{\cal B}^2.\tag{4}$$
It is easy to check that
$${\cal B}_i^a\to \sum_{b=1}^2 M^a{}_b {\cal B}_i^b\tag{5}$$
implements the sought-for $SO(2)$ transformation on the magnetic and electric fields $\vec{B}$ and $\vec{E}$. The two scalar-valued Maxwell equations
$$\vec{\nabla}\cdot\vec{B}~=~0
\qquad\mathrm{and}\qquad
\vec{\nabla}\cdot\vec{E}~=~0 \tag{6}$$
are identically satisfied, because
$$\vec{\nabla}\cdot\vec{\cal B}^a
~=~\vec{\nabla}\cdot(\vec{\nabla}\times\vec{\cal A}^a)
~=~0, \qquad a=1,2.\tag{7}$$
II) Lagrangian density. The Lagrangian density ${\cal L}$ is [2]
$$
{\cal L} ~=~ \frac{1}{2}\sum_{a,b=1}^2 \vec{\cal B}^a \cdot\left( \epsilon_{ab} \frac{\partial \vec{\cal A}^b }{\partial t} - g_{ab} \vec{\cal B}^b\right), \tag{8}
$$
where $\epsilon_{ab} = - \epsilon_{ba}$ is the Levi-Civita tensor in 2 dimensions (with $\epsilon_{12} = 1$). It is easy to check that the Lagrangian density ${\cal L}$ is manifestly invariant under global $SO(2)$ transformations, because both $\epsilon_{ab}$ and $g_{ab}\equiv \delta_{ab}$ are invariant tensors for $SO(2)$. The action is by definition
$$ S[{\cal A}]~=~\int d^4x \ {\cal L}. \tag{9}$$
Extremizing the action with respect to $2\times 3=6$ gauge potential fields ${\cal A}^a_i$ produces $6$ Euler-Lagrange equations. In detail, an arbitrary infinitesimal variation ${\cal A}^a_i\to {\cal A}^a_i+\delta{\cal A}^a_i$ induces a change $\delta{\cal L}$ in the Lagrangian density
$$\begin{align}
\delta{\cal L} ~\sim~& \sum_{a,b=1}^2\delta \vec{\cal B}^a \cdot\left( \epsilon_{ab} \frac{\partial \vec{\cal A}^b }{\partial t} - g_{ab} \vec{\cal B}^b\right) \cr
~\sim~& \sum_{a,b=1}^2\delta \vec{\cal A}^a\cdot \left( \epsilon_{ab} \frac{\partial \vec{\cal B}^b }{\partial t} - g_{ab} \vec{\nabla} \times \vec{\cal B}^b\right),\end{align}\tag{10}
$$
where the "$\sim$" sign means equality modulo divergence terms. The variation $\delta S=0$ of the action vanishes iff the last parenthesis is zero. This yield precisely the two vector-valued Maxwell equations
$$ \frac{\partial \vec{B}}{\partial t} + \vec{\nabla} \times \vec{E}~=~\vec{0}
\qquad\mathrm{and}\qquad
\frac{\partial \vec{E}}{\partial t} ~=~ \vec{\nabla} \times \vec{B}.\tag{11}
$$
References:
S. Deser & C. Teitelboim, Duality Transformations Of Abelian And Nonabelian Gauge Fields, Phys.Rev.D 13 (1976) 1592.
C. Bunster & M. Henneaux, Can (Electric-Magnetic) Duality Be Gauged?, Phys.Rev.D83 (2011) 045031, arXiv:1011.5889; eq. (2.5).
S. Deser, No local Maxwell duality invariance, Class.Quant.Grav.28 (2011) 085009, arXiv:1012.5109.
J.H. Schwarz & A. Sen, Duality Symmetric Actions, Nucl.Phys.B411 (1994) 35, arXiv:hep-th/9304154.
P. Pasti, D. Sorokin & M. Tonin, Note on manifest Lorentz and general coordinate invariance in duality symmetric models, Phys.Lett. B352 (1995) 59, arXiv:hep-th/9503182.
