I am very new to Time Series Analysis (together with R). I have been practising with some very simple datasets to understand how to decompose the time series into Trend, Seasonal and Error components and then check whether the Error component is Gaussian or not. Till this point, everything is pretty much clear to me. However, after the decomposition, how do I fit the trend and seasonal component (using R)? I am confused on 'to what' do I fit my trend and seasonal component? I mean what should my predictors be? For the trend, I can use the time (specifically for the example below, the four quarters) perhaps? And as such, what does the fitting of seasonal component imply and how can it be achieved?

Below is one of time series datasets that I used. This time series is called 'jj' and is present in the 'astsa' package. The frequency of this time series is quarterly.:

        > jj
             Qtr1      Qtr2      Qtr3      Qtr4
      1960  0.710000  0.630000  0.850000  0.440000
      1961  0.610000  0.690000  0.920000  0.550000
      1962  0.720000  0.770000  0.920000  0.600000
      1963  0.830000  0.800000  1.000000  0.770000
      1964  0.920000  1.000000  1.240000  1.000000
      1965  1.160000  1.300000  1.450000  1.250000
      1966  1.260000  1.380000  1.860000  1.560000
      1967  1.530000  1.590000  1.830000  1.860000
      1968  1.530000  2.070000  2.340000  2.250000
      1969  2.160000  2.430000  2.700000  2.250000
      1970  2.790000  3.420000  3.690000  3.600000
      1971  3.600000  4.320000  4.320000  4.050000
      1972  4.860000  5.040000  5.040000  4.410000
      1973  5.580000  5.850000  6.570000  5.310000
      1974  6.030000  6.390000  6.930000  5.850000
      1975  6.930000  7.740000  7.830000  6.120000
      1976  7.740000  8.910000  8.280000  6.840000
      1977  9.540000 10.260000  9.540000  8.729999
      1978 11.880000 12.060000 12.150000  8.910000
      1979 14.040000 12.960000 14.850000  9.990000
      1980 16.200000 14.670000 16.020000 11.610000

        > decomp_JJ=stl(log(jj),s.window = 4)
        > plot(decomp_JJ)
[![STL Plot for JJ ][1]][1]
        
 Any help on this is much appreciated!      
        
        


  [1]: https://i.sstatic.net/8zWES.png