How to plot parameters timeseries in rstan

I have an already rstan fitted object named "fit", of the kind:

Inference for Stan model: MySecondModel.
1 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

                                 mean se_mean    sd        2.5%         25%         50%         75%       97.5% n_eff Rhat
step_size                        0.01    0.00  0.00        0.01        0.01        0.01        0.01        0.01   493 1.00
theta_init                       0.18    0.00  0.02        0.14        0.17        0.18        0.19        0.22  1167 1.00
theta_steps_unscaled[1]         -0.02    0.02  0.91       -1.81       -0.60        0.00        0.58        1.73  1453 1.00
theta_steps_unscaled[2]         -0.02    0.03  0.94       -1.95       -0.65        0.00        0.57        1.81  1382 1.00
theta_steps_unscaled[3]          0.14    0.03  0.95       -1.73       -0.51        0.15        0.80        2.01  1224 1.00
theta_steps_unscaled[4]         -0.34    0.02  0.91       -2.05       -0.95       -0.36        0.26        1.44  1499 1.00
theta_steps_unscaled[5]         -0.32    0.02  0.88       -2.02       -0.94       -0.31        0.31        1.32  2079 1.00
theta_steps_unscaled[6]         -0.99    0.02  0.86       -2.69       -1.56       -0.98       -0.41        0.66  1759 1.00
theta_steps_unscaled[7]          0.23    0.03  0.94       -1.59       -0.40        0.24        0.88        2.06  1411 1.00
theta_steps_unscaled[8]         -0.34    0.02  0.89       -2.13       -0.94       -0.38        0.24        1.34  1774 1.00
theta_steps_unscaled[9]         -0.78    0.02  0.89       -2.49       -1.39       -0.75       -0.17        0.86  1509 1.00
theta_steps_unscaled[10]        -0.66    0.02  0.87       -2.34       -1.28       -0.71       -0.05        1.05  1928 1.00

I want to access the whole vector "theta_steps_unscaled", first to eventually debug the chain, and finally plot for example its mean/mode as a timeseries. I have been looking for something like regex_pars but there does not seem to be any for the "stan" objects.

The best result I was able to achieve was the one below , but this is not what I actually wanted.

I have tried also this, to analyze the chain:

traceplot(fit, "theta_steps_unscaled[6]")

Error in mcmc.list(x) : Arguments must be mcmc objects

And finally casting it as a mcmc:

traceplot(as.mcmc(fit, "theta_steps_unscaled[6]"))

While this raises:

Error in xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 'list' object cannot be coerced to type 'double'

Please note my last goal is to plot the whole vector theta_steps_unscaled[1:some_index], obtaining a similar . Namely with Credibility Intervals, mode/mean and other bayesian parameters.

tidybayes is a handy way to pull out samples in a convenient format. (It has convenient plotting functions too; this vignette gives an overview.) I would use tidybayes to extract your samples and ggplot2 to plot them in the way you want.

First, extract the sampled values:

library(tidyverse)
theme_set(theme_bw())
library(tidybayes)

samples.df = fit %>%
  spread_draws(theta_steps_unscaled[t]) %>%
  ungroup()

Next, plot them with the step on the x-axis and the credible intervals on the y-axis:

samples.df %>%
  group_by(t) %>%
  summarise(median = median(theta_steps_unscaled),
            upper.95 = quantile(theta_steps_unscaled, 0.975),
            lower.95 = quantile(theta_steps_unscaled, 0.025)) %>%
  ggplot(aes(x = t, y = median)) +
  geom_line(size = 1, color = "blue") +
  geom_ribbon(aes(ymin = lower.95, ymax = upper.95),
              fill = "blue", alpha = 0.1) +
  scale_x_continuous("step", breaks = 1:10) +
  scale_y_continuous(expression(theta))

The plot will look something like this: