The dynamics and predictability of stratospheric low-frequency variability in the Northern Hemisphere winter are examined using a two-dimensional (2D) phase space spanned by the leading empirical orthogonal functions of the 10-hPa geopotential height field. The 2D phase space represents the variation of the strength of the polar night jet and the amplitude of zonal wavenumber-1 eddy components. A linearized nonstationary fluctuation dissipation relation (NFDR) is developed based on the deterministic drift vector and the stochastic diffusion tensor estimated from a reanalysis dataset. The authors find that the solution of the linearized NFDR with an optimal data sampling time interval for estimating the drift vector and the diffusion tensor provides a good representation in the phase space of the inhomogeneous distribution of the forecast spread of the operational ensemble forecast conducted by the Japan Meteorological Agency. In particular, the linearized NFDR captures the local maximum of the forecast spread during the onset period of the major stratospheric sudden warming events.