HereTo achieve a symmetric start and end transition necessitates using a filter with a symmetric impulse response, as the resulting step response is the integration of the impulse response. A symmetric impulse response is a linear phase filter, thus the only way to truly do this is with an FIR filter. However linear phase filters can be approximated with IIR solutions, and at the end I show a post-processing "trick" in which an IIR filter could be used if post-processing was feasible for the application.
That said, here are some options to consider:
The first plot is the input step.
The second plot is the result with the exponential moving average filter as given by the OP.
The third plot is the result of using two moving average filters in cascade. (A moving average filter is an FIR filter with unity gain coefficients and then scaled by the sum of the coefficients, so passing through this filter twice or convolving the coefficients to be a single filter). You can also achieve a moving average filter by cascading an accumulator with a difference filter (as a CIC) as follows:
I used two moving average filters each over 40 samples for the resulting plot I gave above so the difference if done as a CIC would be $1-z^{-40}$.
The fourth plot is the result of using a Gaussian Filter. This is the result of a FIR filter with coefficients given by a Gaussian. You can also get that by repeatably passing the signal through a two tap FIR filter with unity gain coefficients scaled by 1/2 ([0.5, 0.5]), which will approach a Gaussian as given by the Central Limit Theorem.
Another option is to use a standard window as the coefficients for an FIR filter (any of the windows that are available in Matlab, Octave, or scipy.signal would be suitable creating such transitions, resulting in a symmetric weighted moving average which is what is needed).
Below is a demonstration with the Hann window and Kaiser window. The transition duration is set by the length of the window used. I like Kaiser as you can adjust the taper directly to match that of most of the other windows just by changing $\beta$.
And a final option if it is feasible to post process the waveform is to pass the waveform through the Exponential Moving Average Filter, and then time-reverse the result, and pass it through the same filter again. This results in a non-causal linear phase IIR filter solution, which provides the desired start and end transition. The command filtfilt available in Matlab, Octave and Python's scipy.signal provides this functionality with the result given below for the OP's EMA filter coefficients. Notice how the result is non-causal, and for that reason I shifted the start transition further to the right to capture the full start of the resulting output:
