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.
1) 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\to \sum_{b=1}^2M^a{}_b {\cal A}^b_i $$
in the 2-dimensional fundamental representation of $SO(2)$, where
$$M^a{}_b =\left[\begin{array}{cc} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) &\cos(\theta)\end{array}\right] \in SO(2), \qquad \sum_{b,c=1}^2 (M^{t})_a{}^b g_{bc} M^c{}_d = g_{ad}, \qquad g_{ab}\equiv \delta_{ab}.$$
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, $$
where
$$ \vec{B}\equiv\vec{\cal B}^1 \qquad \mathrm{and} \qquad \vec{E}\equiv \vec{\cal B}^2.$$
It is easy to check that
$${\cal B}_i^a\to \sum_{b=1}^2 M^a{}_b {\cal B}_i^b$$
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 $$
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.$$
2) Lagrangian density. The Lagrangian density ${\cal L}$ is
$$
{\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),
$$
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}. $$
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
$$
\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)
~\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),
$$
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}.
$$
References:
1) S. Deser and C. Teitelboim, "Duality Transformations Of Abelian And Nonabelian Gauge Fields", Phys.Rev.D 13 (1976) 1592.
2) C. Bunster and M. Henneaux, "Can (Electric-Magnetic) Duality Be Gauged?", Phys.Rev.D83 (2011) 045031, arXiv:1011.5889.
3) S. Deser, "No local Maxwell duality invariance", Class.Quant.Grav.28 (2011) 085009, arXiv:1012.5109.