Why does this reconstruction produce a phase shift from the original signal?

Question

After an FFT of a signal is done, it is plotted as in the image below, with the original signal is on the first subplot.

Using the magnitude and phase data of each frequency, it's reconstructed, and the produced signal is as the image below.

The reconstructed is signal is good except it's shifted 4 milliseconds earlier from the the original signal. I also has done this with another signal, but this 4 milliseconds shift is there too. Why does this happen?

Additional question: I used this for harmonic analysis of 50 Hz fundamental frequency (electric power system frequency). I had a slightly hard time determining which frequency of each harmonic (50*n Hz) frequency to use because, from the FFT result, the phase angle is θ at (50*n∓0.001) Hz but suddenly becomes θ±180° at (50*n±0.001) Hz. For example, at 149.999Hz the phase angle is 223.3° but at 150.001Hz the phase angle is suddenly 403.3°. Is this hard time determining which angle of each frequency to use for reconstruction, normal?

By the way, this is the Matlab code I use to do the FFT and the plotting.

Fs = 1/(8e-6);                % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2^21;                     % Length of signal
t = (0:L-1)*T;                % Time vector

SampleNoBallast;              % The original signal

NFFT = 2^nextpow2(L);         % Next power of 2 from length of y
Y = fft(y,NFFT)/L;            % The Fast Fourier Transform producing FFT complex
f = Fs/2*linspace(0,1,NFFT/2+1);     % Frequencies to plot
P = rad2deg(unwrap(angle(Y)));       % Phase degrees from FFT complex

% Plotting the original signal
subplot(3,1,1); plot(t(1:12500),y(1:12500))
title('Arus Masukan LED T8 Opple 18 W tanpa Ballast')
ylabel('arus (mA)')
xlabel('waktu (s)')
grid on
set(gca,'ButtonDownFcn','selectmoveresize');
% Plot single-sided amplitude spectrum.

% Plotting the magnitude of each frequency
subplot(3,1,2); plot(f(1:17000),2*abs(Y(1:17000))) 
title('Hasil FFT')
xlabel('Frekuensi (Hz)')
ylabel('Amplitudo (mA)')
set(gca,'ButtonDownFcn','selectmoveresize');

% Plotting the phase of each frequency    
subplot(3,1,3); plot(f(1:17000),P(1:17000)) 
title('Sudut Komponen Harmonik')
xlabel('Frekuensi (Hz)')
ylabel('Sudut (derajat)')
set(gca,'ButtonDownFcn','selectmoveresize');

Answer 1

The shift is due to using an FFT with a different length than the length of the data, and likely using a non-symmetric arrangement of zero-padding to increase that original length to the zero-padded length.

Non-symmetric zero-padding rotates the phase results of an FFT, spiraling across result bins.