Does the sun have a vertical movment? As the Sun moves along with the planets does the sun exhibit vertical movement? That is, does it traverse in different planes?
Does the sun have a vertical movment?


Sure it does!

People who follow solar system bodies with the highest accuracy, including NASA's JPL, use the solar system barycenter (center of mass) for the origin of the plane of the ecliptic.

Since the big planets (Jupiter et al.) have some inclination, they pull the Sun up and down out of the plane of the Ecliptic by tens of thousands of kilometers. It's a small fraction of it's 600,000 kilometer radius, but it certainly happens.

You can see the data from the JPL Horizons website, or if you like Python, you can get virtually the same numbers by running the Skyfield package.

From Horizions I downloaded the position of the Sun and the four major planets of the solar system in 10 day steps from the years 1900 to 2100. Below I plot the three components X, Y and Z in the J2000.0 Solar system barycentric coordinates using the ecliptic as the reference frame. The black line is the motion of the Sun, and the color lines are for the four planets, but I have multiplied each by the the mass of the planet divided by the Mass of the Sun.

In the second plot, I added the four together, and you can see that together they mirror the motion of the Sun. 

The last two show the movement of the Sun in 3D (two views). Note that the vertical scale is much smaller than the horizontal scale. 



Horizons setup:





below: Black thick line: motion of the Sun in barycentric coordinates in the J2000.0 X, Y and Z directions, relative to the ecliptic. Four thinner color lines: motions of the four largest planets, each weighted by its mass relative to the Sun's mass. Short period lines are Jupiter and Saturn with Jupiter the largest amplitude, long periods are Uranus and Neptune, with Neptune the largest amplitude.



below: Black thick line: motion of the Sun in barycentric coordinates in the J2000.0 X, Y and Z directions, relative to the ecliptic. Thin red line: sum of the motions of the four largest planets, each weighted by its mass relative to the Sun's mass.



below: Two views  of the 3D motion of the Sun in barycentric coordinates between the years 1900 and 2100. Units are kilometers, and note that the scales are different for each axis.





class Body(object):
    def __init__(self, name):
        self.name = name

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

names = ('Sun', 'Jupiter', 'Saturn', 'Uranus', 'Neptune')

bodies = []
for name in names:
    fname = 'sun motion ' + name + ' horizons_results.txt'
    body  = Body(name)
    bodies.append(body)

    with open(fname, "r") as infile:
        lines = infile.read().splitlines() 

    iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
    iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

    print (name, iSOE, iEOE, lines[iSOE], lines[iEOE])

    lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
    JD  = np.array([float(line[0]) for line in lines])
    pos = np.array([[float(item) for item in line[2:5]] for line in lines])
    body.lines = lines
    body.JD    = JD
    body.pos   = pos.T.copy()

years = 2000 + (JD - JD[3652])/365.2564
Sun, Jupiter, Saturn, Uranus, Neptune = bodies

for body in bodies:
    body.rsun = np.sqrt(((body.pos - Sun.pos)**2).sum(axis=0))

planets = [b for b in bodies if 'sun' not in b.name.lower()]

massdict = {'Jupiter':1.898E+27, 'Saturn':5.683E+26, 'Uranus':8.681E+25,
            'Neptune':1.024E+26, 'Sun':1.989E+30}

for b in bodies:
    b.mass = massdict[b.name]

if True:
    plt.figure()
    for i, thing in enumerate(Sun.pos):
        plt.subplot(3, 1, i+1)
        plt.plot(years, thing, '-k', linewidth=1.5)
    for p in planets:
        for i, thing in enumerate(p.pos):
            plt.subplot(3, 1, i+1)
            thing = thing * (p.mass/Sun.mass)
            plt.plot(years, thing)
    plt.show()

all_four = Body('all_four')
all_four.pos = np.array([p.pos*(p.mass/Sun.mass) for p in planets]).sum(axis=0)

if True:
    plt.figure()
    for i, thing in enumerate(Sun.pos):
        plt.subplot(3, 1, i+1)
        plt.plot(years, thing, '-k', linewidth=1.5)
    for i, thing in enumerate(all_four.pos):
        plt.subplot(3, 1, i+1)
        plt.plot(years, thing, '-r')
    plt.show()

if True:    
    fig = plt.figure(figsize=[12, 10])  # [12, 10]
    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = Sun.pos
    ax.plot(x, y, z)

    plt.show()