Simulating Phase Noise: How to do it?

I have a system that I am simulating, and I want to add phase noise to the process. I didn't know how to do it, so I have been conducting an literature review.

I found a good reference with an "easy to understnad" example. This example is on page 145 and 146, where there are a pair of equations:

$$osc(t[n])=\sin(\omega_0 t[n] + \Delta \phi(t[n])$$ $$\Delta \phi(t[n])=\phi_{synch}(t[n])+\sum^n_{k=1}\phi_{acc}(t[k])$$

In these equations, $osc(t[n])$ is the oscillator output, $\omega_0$ is the frequency of oscillation, $\Delta \phi(t[n])$ is the phase noise, and $\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$ represent the synchronous and accumulating components of the phase deviations.

The authors used these equations in conjunction with some matlab code to publish a nice graph on page 146:

This is exactly what I need! However, I am confused by the matlab code, and by the exact nature of $\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$. Specifically, the authors say:

$\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$ are random variables with Gaussian distributions.

These are my questions:

And what purpose do the filters in the code (shown below) serve?

Bonus quote from this fantastic book:

In order to derive an expression for $\Delta \phi(t[n])$, we note that an oscillator is an autonomous system in which phase deviations accumulate over time. However, the inherent accumulation of phase deviations does not affect the noise contributions of external sources that are added directly to the output of the oscillator , as shown in [4]. Hence ,the phase of the output signal of an oscillator can include a noise component which is not accumulated.This noise component can be considered synchronous with the output of the oscillator. Accumulating and synchronous phase deviations are related to accumulating and synchronous jitter...

Code from appendix A.4 page 173:

clc
close all
clear all

tstep = 1;
nsamples = 2^23
time = 0:tstep:nsamples -1;
period = 2^9;
ThreeSigma = period/20;
Sigma = ThreeSigma /3;
varjitter = Sigma^2;
Nperiods = length(time)/period +10;

jitterWhite = Sigma*randn(1,nsamples);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% FILTER DESIGN: Single POLE %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms z z1 outz inz p1;
fs = tstep/period;
s = 2 * fs * (z-1)/(z+1);
Hz = p1/(s+p1);
[N,D] = numden(Hz);
Cnum = coeffs(N, z);
Cden = coeffs(D, z);
Cnum1 = Cnum(1);
Cnum2 = Cnum(2);
Cden1 = Cden(1);
Cden2 = Cden(2);
fsval = 2*pi*fs; % sampling frequency in the discrete time domain
p1val = fsval /100; % frequency of the pole
Cnum1valextra = double(subs(Cnum1 , {p1, fs},{p1val , fsval}));
Cnum2valextra = double(subs(Cnum2 , {p1, fs},{p1val , fsval}));
Cden1valextra = double(subs(Cden1 , {p1, fs},{p1val , fsval}));
Cden2valextra = double(subs(Cden2 , {p1, fs},{p1val , fsval}));
Aextra = (Cnum2valextra/ Cden2valextra);
Bextra = (Cnum1valextra/ Cden2valextra);
Cextra = - (Cden1valextra/ Cden2valextra);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FILTER: Implementation %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jitterFiltered = zeros(1,length(jitterWhite));
for n = 2:1: length(jitterWhite)
jitterFiltered(n) = (Cnum2valextra/ Cden2valextra) * jitterWhite(n) + (Cnum1valextra/ Cden2valextra) * jitterWhite(n-1) - (Cden1valextra/ Cden2valextra) * jitterFiltered(n-1) ;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% SECOND FILTER: Design and Implementation %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jitterWhiteNew = Sigma*randn(1,nsamples);
syms z z1 outz inz p1;
fs = tstep/period;
s = 2 * fs * (z-1)/(z+1);
Hz = p1/(s+p1);
[N,D] = numden(Hz);
Cnum = coeffs(N, z);
Cden = coeffs(D, z);
Cnum1 = Cnum(1);
Cnum2 = Cnum(2);
Cden1 = Cden(1);
Cden2 = Cden(2);
fsval = 2*pi*fs; % sampling frequency in the discrete time domain
p1val = fsval /10000; % frequency of the pole (Hz)
Cnum1valextra = double(subs(Cnum1 , {p1, fs},{p1val , fsval}));
Cnum2valextra = double(subs(Cnum2 , {p1, fs},{p1val , fsval}));
Cden1valextra = double(subs(Cden1 , {p1, fs},{p1val , fsval}));
Cden2valextra = double(subs(Cden2 , {p1, fs},{p1val , fsval}));
Aextra = (Cnum2valextra/ Cden2valextra);
Bextra = (Cnum1valextra/ Cden2valextra);
Cextra = - (Cden1valextra/ Cden2valextra);
jitterFilteredNew = zeros(1,length(jitterWhite));
for n = 2:1: length(jitterWhite)
jitterFilteredNew(n) = (Cnum2valextra/ Cden2valextra) * jitterWhiteNew(n) + (Cnum1valextra/Cden2valextra) * jitterWhiteNew(n-1) - (Cden1valextra/ Cden2valextra) * jitterFilteredNew(n-1) ;
end
jitter = 0.01*( 0.01* cumsum(jitterFilteredNew) + 0.01* jitterFiltered );
sinsignalNoisy = sin(2*pi*(1./ period).*time+ jitter);